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Preface 



One of the key features of QCD is asymptotic freedom. While the theory is amenable 
to perturbation theory at large momenta, it is non-perturbative for energy scales < 1 
GeV and lattice QCD is the only known method for first principles calculations in 
this regime. The running of the coupling with momentum scale immediately implies 
the existence of different states of nuclear matter at asymptotically high densities or 
temperatures: when the coupling on the scales of temperature T or baryon chemical 
potential /is is sufficiently weak, confinement gets lost. At a critical temperature 
200 MeV, QCD predicts a transition between the familiar confined hadron physics and 
a dcconfined phase of quark gluon plasma (QGP). At the same temperature, chiral 
symmetry gets restored. A thermal environment with sufficiently high temperatures 
for a QCD plasma has certainly existed during the early stages of the universe, which 
passed through the quark hadron transition on its way to its present state. On the 
other hand, for high densities and low temperatures, there is an attractive interaction 
between quarks to form Cooper pairs, and exotic non-hadronic phases such as a colour 
superconductor have been predicted. Such physics might be realised in the cores of 
compact stars. 

Current and future heavy ion collision experiments are attempting to create hot 
and dense quark gluon plasma at RHIC (BNL), LHC (CERN) and FAIR (GSI). These 
studies will have a bearing far beyond QCD in the context of early universe and astro- 
particle physics. Many other prominent features of the observable universe, such as 
the baryon asymmetry or the seeding for structure formation, have been determined 
primordially in hot plasmas described by non-abelian gauge theories. The QCD plasma 
serves as a prototype also for those, since it is the only one we can hope to produce 
in laboratory experiments. There is thus ample motivation to provide controlled the- 
oretical predictions for the physics of hot and dense QCD. As we shall see, this turns 
out to be even more challenging than QCD in the vacuum. 

In these lectures, the focus is on basic concepts and methods, with a few exemplary 
results for illustration only and a very incomplete list of references. For summaries of 
the latest results and literature lists, see the annual proceedings of the LATTTICE 
conferences. 
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Aspects of finite temperature field 
theory in the continuum 



1.1 Statistical mechanics reminder 

We wish to describe a system of particles in some volume V which is in thermal 
contact with a heat bath at temperature T. Associated with the particles may be a set 
of conserved charges A'^^, i = 1, 2, . . . (such as particle number, electric charge, baryon 
number etc.). In statistical mechanics there is a choice of ensembles to describe this 
situation. In the canonical ensemble, V and the Ni are kept fixed while the system 
exchanges energy with the heat bath. In the grand canonical description, also exchange 
of particles with the heat bath is allowed. Of course, both ensembles provide the same 
description for physical observables in the thermodynamic limit, — > oo, but for 
specific situations one or the other may be more appropriate. In quantum field theory, 
the most direct description is in terms of the grand canonical ensemble. Its density 
operator and partition function are given as 

p = c-iiH-,,N,) ^ Z = Trp, Tr(...) = ^(n|(...)|n) , (1.1) 

n 

where /i^ are chemical potentials for the conserved charges, and the quantum mechan- 
ical trace is a sum over all energy eigenstates of the Hamiltonian. Thermodynamic 
averages for an observable O are then obtained as (O) = Z^^'Tr{pO) . 

From the partition function, all other thermodynamic equilibrium quantities follow 
by taking appropriate derivatives. In particular, the (Helmholtz) free energy, pressure, 
entropy, mean values of charges and energy are obtained as 



F=-TlnZ, d{T\nZ) 
^^a(TlnZ) = 

9y ' E = -pV + TS + p,N,. (1.2) 

_ d{T\nZ) 

dT ' 

Since the free energy is known to be an extensive quantity, F = fV, and we are 
interested in the thermodynamic limit, it is often more convenient to consider the 
corresponding densities, 

f=y, P = -L s = -, n, = -, 6=-. (1.3) 
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1.2 QCD at finite temperature and quark density 

Let us now consider the grand canonical partition function of QCD. The derivation of 
its path in tegral representation is discussed in detail in the textbooks ( Kapusta and 
Gale, 2006rald requires discretised Euclidean time followed by a continuum limit. We 
shall thus derive it more conveniently in its lattice version in Sec. 12.11 and just quote 
the result here, 

Z{V,^Jif,T■g,mf)=TT(c-'^"-^'f'^f^'^^ = J DA D^} e~^^^^-^ e-'^f^^-^-^-^ 

(1-4) 

with the Euclidean gauge and fermion actions 



l/T 

j dxo j ^Tr F^i,(a;)Fpi,(a;), 



V 

Sf[%i),'ii),Af,]^ j dxo (fx ^i)f{x) {^f,Df, + mj - fif^o)iJf{x). (1.5) 

V 

The index / labels the different quark flavours, and the covariant derivative contains 
the gauge coupling g, 

D^ = {d^^igA^), A^=T-Al{x), a^l,...N^-l, 

F^,{x) = ^[Df,,D,] . (1.6) 

The thermodynamic limit is obtained by sending the spatial three- volume V — )• oo. 
The difference to the Euclidean path integral at T = is that the temporal direction 
is compactified and kept finite, i.e. the theory lives on a torus whose compactification 
radius defines the inverse temperature, l/T. The path integral is to be evaluated with 
periodic and anti-periodic boundary conditions in the temporal direction for bosons 
and fermions, respectively, 

A^(t,x) - A^(T + i x), V(r,x) = -V(r + ^,x) , (1.7) 

which ensures Bose/Einstein statistics for bosons and the Pauli principle for fermions. 
Clearly, the path integral for vacuum QCD on infinite four-volume is smoothly recov- 
ered from this expression for T — > 0. 

The partition function depends on the external macroscopic parameters T, y,/x/, 
as well as on the microscopic parameters like quark masses and the coupling constant. 
The conserved quark numbers corresponding to the chemical potentials /x/ are 

Qf = i'flo^f ■ (1-8) 



Since the QCD phase transition happens on a scale 200 MeV, we neglect the c, b, t 
quarks or treat them non-relativistically when needed. We will thus consider mostly 
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two and three flavours of quarks, and always take m„ = m^. The case nis = m.u^d is 
then denoted by Nf ~ 3, while Nf — 2 + 1 implies rrig ^ rriu^d- Furthermore, we will 
couple all flavours to the same chemical potential fi unless otherwise stated. 

The only, but fundamental, difference compared to the T = path integral is due 
to the compactness of the time direction. For example, the Fourier expansion of the 
fields on a finite volume = is 

oo 

^ n— — OO p 

^ oo 

ij{r,^)^— ^n = {2n+l)7TT . (1.9) 

^ n— — oo p 

The allowed momenta are pi = {2Trni)/L, where rii G 1. The normalisation factors 
in front are chosen such that the Fourier modes, called the Matsubara modes, are 
dimensionless. In the thermodynamic limit the momenta become continuous 



(1.10) 



V ^ J (27r)3' 

ni,n2,ri3 ' 

but the Matsubara frequencies w„ stay discrete due to the compactness of the time 
direction. Note that bosons/fermions have even/odd Matsubara frequencies, respec- 
tively, to ensure the bondary conditions Eq. (|1.7p . We thus obtain modified Feynman 
rules ( [Kapusta and Gale, 2006| compared to the vacuum. The zero components of 
four vectors contain the Matsubara frequencies w„, and the four- momentum integra- 
tion associated with internal lines gets replaced by a 3d integral and a Matsubara sum. 



E 

n— — oo 



(1.11) 



(27r) 



1.3 Perturbative expansion 

Similar to T = 0, we can take the path integral as starting point for a perturbative 
expansion. For some generic field 0, we proceed just as in the vacuum and split the 
action in a free and an interacting part, S = {Sq + Si), in order to expand in powers 
of the interaction, and hence the coupling constant, 

Z = N I D(j) e-(^«+^') = N I D<p e'^" E ^T"*^* • (l-l^) 

Thus we get for the log of the partition function 

In Z = In Zo+ln = In [n J e"^") +ln (^1 + E ^ J^ofJ-sf ) ' (^-l^) 



We shall evaluate the ideal gas part, InZo, in the next section. The corrections to the 
non-interacting case are the sum of all loop diagrams without external legs. When 
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evaluating loop diagrams, UV divergences are encountered and the renormalisation 
program has to be performed. Whatever regularisation and renormalisation is neces- 
sary and sufficient at zero temperature is also necessary and sufficient at finite tem- 
perature. This is not surprising, since the ultraviolet structure of the theory, i.e. the 
microscopic short distance regime, is unchanged by the introduction of macroscopic 
parameters like temperature and baryon chemical potential. 

On the other hand, the infrared structure of the theory does get changed. This is 
easily seen by considering the inverse propagator for a bosonic degree of freedom. 



The Matsubara frequencies act like effective thermal masses ~ T for all modes with 
n 0. In the case of fermions, there are only non-zero modes. Hence, all fermionic and 
all non-zero bosonic modes are infrared-safe even in the limit of vanishing bare mass, 
m ~ Q. By contrast, for the bosonic zero mode the inverse propagator is -I- m^, 
which is identical to that of a 3d field theory. Thus, 4d Yang-Mills theory at finite 
temperature contains in the zero mode sector the 3d Yang-Mills theory, which is a 
confining theory. Its propagator is infrared divergent for m = 0, and the divergence is 
worse than in 4d. This points at non-perturbative behaviour and is at the heart of the 
Linde problem of finite temperature perturbation theory. Sec. 11.61 

Let us briefly discuss qualitatively the self-energy corrections to the gluon propa- 
gator. When computing the self energy diagrams, one finds that for the colour electric 
field ^0 a gluon mass is generated, the electric or Debye mass. To leading order it is 



Thus, colour electric fields get screened by the medium at finite temperature, whereas 
for colour magnetic fields Ai one finds the corresponding magnetic mass rriM = at 
this order. However, there are contributions to ttim ~ g^T starting at the two-loop 
level. As we shall see in Sec. II. 6[ all loops contribute equally to the coefficient, such 
that the calculation of the magnetic screening mass is an entirely non-perturbative 
problem. 

1.4 Ideal gases 

Ideal, i.e. non-interacting, gases of particles are important model systems to guide our 
intuition. It is therefore instructive to see how their thermal properties are derived from 
the path integral. Let us consider a bosonic system of real scalar fields, for simplicity. 
After Fourier transformation, the action reads 



(1.14) 




(1.15) 




(1.16) 



where we abbreviate uj ~ -^/p^-l-m^ and (f>n{p) — 4'-n{p) for a real scalar field. The 
partition function then factorises into a product over all Matsubara modes, 
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oo 

Zo = N U [] / d0„ cxph2y2K'+^')'^"(P)Cb)] 

n— — oo p 

=ivn 
xn 

n>0'' 



(1.17) 



n>0 ' 



iV'n(2vr)i 



/2 f ^ 



y2 



Tl>0 



4 T2 



" n n 



n— — oo p 



T2 



where we have changed to polar coordinates in the third hne and absorbed constant 
numerical factors into the normalisation in front. Thus, just as for zero temperature, 
the partition function of the free theroy can be formally written as 



Za = N D(t> e-^(*) = iV"(det ^-^y^'^ 



(1.18) 



where = (a;^ + uj'^)/T'^ is the inverse propagator in momentum space. 

Since N" is (V, r)-independent, it will not contribute to thermodynamics, Eqs. (|1.2p . 
and may be dropped. (It will contribute e.g. to the entropy as an additive constant, 
which we are allowed to set to zero by the third law of thermodynamics). Thus we 
obtain 



InZo 



n— — oo p 

The Matsubara sum is performed using the formulae 



^■2 



(1.19) 



In 



UJ 



E 



1 



n2 + (^)2 



612 + (27rn)2 
2 



+ ln(l + (27m)2) 



1 



(1.20) 



leading to the expression 
In Zo = 



E 



V 



(2^ 



de 



2T 



2 e»-l 



ln(l -e-T^ 



T-indcp. 



(1.21) 



One observes that the integral over the first term diverges in the UV since to ~ |p|. 
This however is familiar from the quantum mechanical harmonic oscillator. The term 
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InZo = -V^ / In(l-c-T) . (1.23) 



corresponds to the zero point energy and gives a divergent vacuum contribution to the 
energy and pressure for T — > 0. The renormahsation condition is that the vacuum has 
zero pressure, 

Pphys(r)=p(T)-p(T = 0). (1.22) 

The final result then is the familar form of the partition function for a free gas of 
spinless bosons, 

For the massless case, m = 0, the momentum integral can be done exactly and one finds 
the famous result for the pressure of one bosonic degree of freedom at zero chemical 
potential. 

Doing the same calculation for non-abelian gauge fields A^, each field component 
corresponds to one bosonic mode and we have to sum over a = l,...iV^ — 1 and 
/i = 1 . . .4 in Eq. ([LTe]) . Hence we get a factor of 4(7V2 - 1) in front of Eq. (fL23|) . 
Going through a similar sequence of steps for a free Dirac field, one finds instead 



In Zo = 2V 



(2^ 



In 1 + e ~ + In 1 



(1.25) 



The factor of two in front accounts for the two spin states of a fermion, and there 
are two terms from fermion and anti-fermion in the brackets. Recall that for gauge 
theories gauge fixing is necessary in order to invert the two-point function. Thus, 
for the free photon or Yang-Mills gas, two of the four bosonic Lorentz degrees of 
freedom get cancelled by the corresponding ghost contributions, such that wc obtain 
two polarisation states per massless vector particle, as it should be. 
All of this can be summarised by the one-particle partition function 

\nZl{V,T) = r^Vv, J ^ ln(l + 77e-(-'-'^')/^) , (1.26) 

where 77 = —1 for bosons and t] ~ 1 for fermions and ui gives the number of degrees 
of freedom for the particle i. The momentum integration for one massless fermionic 
degree of freedom gives the pressure 

p=-—T\ (1.27) 
^ 8 90 ^ ' 

To compute the pressure of an ideal gas of gluons and massless quarks, we now simply 
have to count the degrees of freedom to obtain (two polarisation states for each of 
{N^ — 1) gluons, two polarisation states per quark, three colours per quark and the 
same for anti-quarks) 

|r=(2(iV2-l)+4iV^/0^- (1-28) 

This is the Stefan-Boltzmann limit of the QCD pressure which is valid for vanishing 
coupling, i.e. in the limit T — > 00. 
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Fig. 1.1 I + 1-loop Fcynman diagram contributing to the pressure. 

1.5 The hadron resonance gas model 

There is another ideal gas system that is useful to model the low temperature behaviour 
of QCD. The modelling assumption is that, at low temperatures when we still have 
confinement, the QCD partition function is close to that of a free gas of hadrons. While 
strong coupling effects are responsible for confinement of the quarks and gluons, the 
interactions between the hadrons are considerably weaker and may even be neglected 
if the gas is sufficiently dilute. In this case, the QCD partition function factorises into 
one-particle partition functions Zl[V,T), 

\nZ{V,T)^Y.^nZ}{V,T) . (1.29) 

i 

The particles to be inserted are the known hadrons and hadron resonances (for QCD 
purposes electroweak decays are to be neglected). Taking the Vi and the energies 
from the particle data booklet, one can supply this formula with hundreds of hadron 
resonances, do the momentum integral and obtain a thermodynamic pressure that 
compares remark ably well with Monte Carlo simulations of QCD ( Karsch, Redlich 
and Tawfik, 2003'). We shall see in Sec. 12.91 that this is no accident, but that the 
hadron resonance gas model can actually be derived from lattice QCD as an effective 
theory for the strong coupling regime. 

1.6 The Linde problem 

As an example to illustrate the Linde problem, consider the I + 1-loop contribution 
to the QCD pressure. Fig. 11.11 with 21 three gluon vertices and 3Z propagators. The 
Matsubara sums over the internal lines, Eq. (|1.11|) . contain a term coming exclusively 
from the zero modes. Dispensing with the index structures, its contribution is given 
by the 3d loop integral 

I^g^' [t I d^p^'p^'iP^+mT''' , (1-30) 

where we have introduced a mass m by hand as an infrared regulator. The momentum 
integrals have to be performed up to a scale T, which appears as an effective UV cut- 
off after doing the Matsubara sum over the other modes. Parametrically, the integral 
then evaluates to 

'g^THn{T/m) for ^ = 3 

g^T\g^T/my-^ for l>3. ^^"^^^ 
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It is therefore infrared divergent for m — >■ 0, Z > 2 and the usual bare perturbation 
theory breaks down. However, as we discussed before, mass scales are generated dy- 
namically by loop corrections to the propagators. Evaluating the integral with those 
mass scales corresponds to an effective resummation of the perturbative series (the 
diagram now contains loop insertions that would make it formally higher order in bare 
perturbation theory). If our internal lines contain Aq fields, we need to do so with 
rriE ^ gT and observe that as a consequence our series contains odd powers of g in 
addition to the logarithms. Summing up mass corrections for the Ai amounts to an 
insertion of rriM ~ g^T, and we see that the gauge coupling drops out of the effec- 
tive expansion parameter entirely for I > 3. Hence, all loop orders contribute to the 
pressure at order g^, which thus is a fully non-perturbative problem. 

The same problem occurs when calculating other quantities, with the order of the 
breakdown depending on the observable. For example, for the electric mass m^; it 
appears already at NLO and for the magnetic mass ttim even at the leading non- 
zero order. Thus at finite temperature, perturbation theory only works to some finite 
order which depends on the observable. Note that this is true no matter how weak 
the coupling g. Even electroweak physics at finite temperatures is inherently non- 
perturbative, and perturbative answers are only useful to the extent that the calculable 
orders are sufficient for a good approximation of physical results. 



2 

The lattice formulation for zero 
baryon density 



2.1 Action and partition function 

Let us now consider the lattice formulation of SU{N) pure gauge theory on a hypcr- 
cubic lattice, x Nr, with lattice spacing a and the standard Wilson gauge action 



Sa[U]^Y. E /? f 1 - ^RcTrt/^y (2.1) 



X l<^<iy<4 

where Up = U^{x)U„{x+afl)Uj^{x+aO)Ul{x) is the elementary plaquette, and the bare 
lattice and continuum gauge couplings are related by /3 = 2N/g'^. As usual, we impose 
periodic boundary conditions in all directions, Ufi{T,x) = C/^(r + iVi-, x), C/^(r, x) = 
(t, X + Ns). The connection between zero and finite temperature physics is most 
easily exhibited by the transfer matrix, which relates the path integral representation 
of a Euclidean lattice field theory to the Hamiltonian formulation. A transfer matrix 
element between two times slices r and r + 1 is given by ( [Montvay and Miinster, 1994[ ), 

T[U,{T + l),U,iT)]=c-''" = J DUoir) exp -L[C/,(t + 1), f/o(T), (/.(r)], (2.2) 

where the action is written as a sum over time slices, 



5g = 5]L[[/,(T + l),[/o(T),[/,(T)], 



L[[/.(T + l),;7o(T),t/,(r)] = ^Li[U,iT + l)] + ^L,mT)] 

+L2[U,{T + l),Uo{T),UdT)] 

Li[C/,(r)] =-^5^ReTrC/p, 

P(t) 

L2[C/.(r + l),;7o(T),t/,(r)] = --^ ^ ReTrUp, (2.3) 

p{t.t + 1) 

and p(T),Ui{T) denote all spatial plaqucttes and links contained in the timeslice r 
with arguments x suppressed, while p(r + 1, r), Uo{t) are timelike plaqucttes and links 
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connecting the timeslices r and r + 1. The partition function is now conveniently 
expressed as 

Z= /n(^C^^(^,x)T[[/,(r + l),CX,(T)])-Tr(T^O-Tr(G-^-'^«). (2.4) 

Note that the periodic boundary condition in the temporal direction is necessary for 
the trace operation in order to have identical states |n) on the time slices 1 and A^,-. 
In this discretised form we can immediately see that Z is equivalent to the partition 
function of a thermal system if we identify 

1 ^ aiV, . (2.5) 

The thermal expectation value of an observable is then 

(O) = Z-iTr(e-f O) = Z'^ ^(n|T^^O|n) = ^"i^^'^'l^X^""" • (2.6) 

As in the continuum, we are interested in the thermodynamic limit and hence — > oo, 
while keeping aNr = 1/T finite. 

In this form we easily see the connection to T = physics: projection on the 
vacuum expectation value is achieved by taking Nt- to infinity, 

= ^1™^ j^^^.aNAE.-E,) ■ (2-7) 

In order to describe our gauge theory at finite temperatures, we simply need to dis- 
pense with this step. In that case the expectation value receives contributions from 
all eigenstates \n) with Boltzmann weights c~^, Eq. p.Op . Hence, all lattices with 
finite Nr (in particular those used for numerical simulations!) correspond to a finite 
temperature T = l/{aNr), and for a description of vacuum physics sufficiently large 
Nr is required. 

On a Euclidean lattice, space and time directions are in principle indistinguishable 
as long as we provide them with equal boundary conditions. For some applications it 
is useful to define a Hamiltonian that translates states in a spatial, say the z-direction, 
defined through a transfer matrix between adjacent z-slices, and write the partition 
function in terms of that Hamiltonian (with U{z) denoting {[/^(z)|/i 7^ 3}), 

r[[/(z + l),C/(z)] =e-"-^% Z = Tr(e-'^^^-^^) . (2.8) 

For thermal physics, we want to take Nx,y,z — > 00 and keep Nr, which is now hidden 
in the definition of the Hamiltonian, finite. It is thus equivalent to the "zero temper- 
ature" {Nz 00) physics of the Hamiltonian H^, which acts on states defined on a 
space with two infinite and one finite, compactified direction. Clearly, has reduced 
symmetry compared to H. In either description, from a calculational point of view fi- 
nite temperature physics on a Euclidean space-time lattice is nothing but a finite size 
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effect: the Boltzmann weighted sums, i.e. thermal effects, become noticeable once the 
temporal lattice size is small enough for the system to be sensitive to the boundary. 

Adding fermions is now the same as for zero temperature and does not change this 
picture. Once a suitable action 



(2.9) 



has been selected, the Gauss integral can be done and, minding the appropriate bound- 
ary conditions in the temporal direction, we end up with the partition function 



Z{Ns,Nr;l3,mf) ^ f DU JJdet M{mf) e" 

C/^(r,x) = C/^(r + iV„x), 
'0(t,x) = -V'(t + Ar^,x) . 



■SAU] 



(2.10) 



For definiteness, the lattice action for Nf degenerate Wilson fermions is given by 



(2.11) 



2.2 Tuning temperature and the continuum limit 



According to Eq. (|2.5p . one way of tuning temperature on the lattice is by choosing 
Nr- But this is not satisfactory as this is only possible in discrete steps, and for real- 
istic lattice spacings these are much too coarse. Hence, common practice in numerical 
simulations is to keep Nr fixed and instead vary the lattice spacing a via the lattice 
coupling, /3 — 2N/g'^{a), thus affecting temperature. This is a marked difference to 
simulations at zero temperature. In particular, simulation points at different temper- 
ature correspond to different lattice spacings and thus have different cut-off effects. 

The relation between a and /3 is given by the renormalisation group. E.g., for the 
Lambda-parameter in lattice regularisation we have to leading order in perturbation 
theory for SU{3) 



aA, = I — 



'bi/2bl 



s_ 

e i2''o 



1 / 2 
bo = ^ 11 - -Nf 



bi 



102 - I 10 + - ) iV/ 



(2.12) 



However, perturbation theory is not convergent for accessible lattice spacings. The way 
out is to express the calculated observables in terms of known physical quantities of 
the same mass dimension. For example, if we want to compute the critical temperature 
of the QCD phase transition, Tc, by keeping Nr fixed and tuning /3, the location of 
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a phase transition will be given as a critical coupling /3c. The critical temperature in 
units of a hadron mass is then 

Tr 1 1 , , 
— — = _ (2.13) 

rriH aciriHNr a{(3c)mHNr 

In order to set the physical scale for Tc, we then have to calculate the zero temperature 
hadron mass in lattice units at the value of the critical coupling, (am^f)(/3c)- 

This procedure is good as long as we are able to simulate physical quark masses. For 
most practitioners, this is not yet the case, and furthermore there are many interesting 
theoretical questions concerning regimes with unphysical quark masses, such as the 
quenched and the chiral limits. In order to set a scale in those cases, one uses quantities 
that display only little sensitivity to the quark mass values. Examples are the string 
tension or the Sommcr scale, 

T 1 ro -ydVir) 

— = ^—, a« 425McV; Tro = r^— = 1.65 . 2.14 
V(T a^aNr aNr dr 

In order to take the continuum limit, we now have to compute expectation values of 
our observables of interest, (0)(/3, m/), in the thermodynamic limit for various lattice 
spacings a. Then we can extrapolate to a — > 0. Since the continuum limit has to be 
taken along lines of constant physics, i.e. keeping temperature and mass ratios fixed, 
this is equivalent to taking — > oo. A continuum extrapolation therefore requires 
simulations on a sequence of lattices with different Nr- For example, at a temperature 
T = 200 McV - 1 fm-\ Eq. with Nr = 4,8,12 implies a w 0.25, 0.125, 0.083 fm. 

2.3 Constraints on lattice simulations and systematic errors 

It is important to realise from the outset that current lattice simulations at finite 
temperature and density are still hampered by sizeable systematic errors and uncer- 
tainties. Let us discuss the origins of those. The Compton wave length of a hadron is 
proportional to its inverse mass ra^ , and the largest of those constitutes the correla- 
tion length of the statistical system. To keep finite size as well as discretisation errors 
small, we need to require 

a<m^^<a7Vs. (2.15) 

For the low T phase we thus need a lattice size of several inverse pion masses. The 
push to do physical quark masses is only just beginning to be a possibility on the most 
powerful machines and with the cheapest actions (i.e. staggered, with Wilson rapidly 
catching up). On the other hand, at high T screening masses scale as mn ~ T, thus 

N-^ 1 NsN-\ (2.16) 

Hence, while we desire large N-r, the spatial lattice size should be significantly larger 
than the temporal one. This limits the directly accesible temperatures to several Tc. 

Another important subject are cut-off effects. Similar to zero temperature, one 
can formally expand lattice observables in powers of the lattice spacing about their 
continuum limit. According to the discussion above, this is tantamount to expanding in 
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the dimcnsionlcss about zero. In order to reduce cut-off effects, one can then apply 
improved actions and operators, again in complete analogy to QCD in the vacuum. 

In order to study thermodynamic behaviour, one is often interested in the tempera- 
ture dependence of certain quantities. This already implies sequences of simulations at 
many /3- values for just one A^,- . If one is interested in a phase transition, one is moreover 
confronted with fluctuations between different phases, increasing correlation lengths 
and critical slowing down, i.e. increased auto-correlation times. For these reasons most 
thermal simulations require orders of magnitude more Monte Carlo trajectories than 
typical zero temperature problems. Simulations are thus necessarily run on compara- 
tively coarse lattices, implying larger systematic errors which are so far less controlled 
than for vacuum physics. 

2.4 The ideal gas on the lattice 

Similar to Sec. 11.41 we can evaluate the ideal gas for a bosonic field on the lattice. 
Starting point is the equation employing the free propagator, 

InZo = - ilndet A = ^TrlnA^^ 
2 2 



where now we have to deal with the lattice momenta, 

;5^=4sin^(^^)+4;f]sin^(^f), 

3 

4ci2 =4;^sin2(^) + (am)^ (2.18) 

Due to the finite lattice spacing the Matsubara sum is only from —Nr/2, ....Nr/2 — 1. 
In order to perform the sum, we employ the formula ( [Kaste and Rothe, 1997| 

^ n=-N^/2 z, 

to the derivative of the sum of logs, 



Nr ^ \ 2 

n 

dL _ I ^ 4 

n 

aiz') 4 



-E- 

A 
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This function has simple poles at z = 2uj'^ + 1 ± 2Cj\/ Qp- + 1. It is now convenient to 
change variables by a) = sinh(ai?/2) to obtain 

dL dL duj^ / 1 

==2 ^rr-T^ - + 1 



d{aE) d{aE) 

L = -|- In (1 - e"^-"^) + 2aE . (2.21) 

We recognise again the vacuum energy contribution, which needs to be subtracted to 
arrive at the final result for the pressure of a bosonic gas on the lattice, 

lnZo = -V |"^^ln(l-e-^-^). (2.22) 

One can now study the approach to the continuum by extending the integration 
range to infinity and expanding the log in small lattice spacing. The leading mistake 
we make by the first step in the positive integration range of the integral is 

^("^ ^ r ^"^^ ~ '"''"''^ ^ (^-^^^ 

and similarly for the negative range. With /(O) ~ and /'(a) cx exp— A^rTi", with 
N-r = l/{aT), this is exponentially small and can be neglected compared to the power 
corrections of the integrand. For the expansion of the integrand, we need the corrections 
to the dispersion relation. Writing 



E{p) = (p) + a^E^^'> (p) + . . . , (p) = Vp^T^, (2.24) 

and expanding both sides of the equation 

sinh2( — ) = ^sm2(^) + 

+ = V f Ml! - + Oia% (2.25) 

4 48 ^l4 48 / 4 ^'^' ^ ' 

j=i ^ / 

one finds for the a^-correction to the dispersion relation 

Changing to dimensionless variables, x = p/T,e ~ E/T, we then have for the expan- 
sion of the pressure 

P ( P\ 2 /" c?^^ e^'^\x) 

^ VT^Jcont J (2^)3 e--*"' W - 1 + • ■ ■ ^^-^^^ 
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This means that the pressure of a free gas of bosons has leading lattice corrections of 
O(a^). For the massless case the integral can again be done and one arrives at 

p 87r2 1 / 1 \ ^ ^ 

= l + ^l^ + T74 > (2-28) 



Pcont 21 V^r 

where Pcont is the continuum result Eq. ()1.24p . 

One can repeat a similar calculation for various fermion actions and discuss their 
cut-off effects. For the free gas in the chiral limit this can be found in ( Hegde, Karsch, 
Laermann and Shcheredin, 2008), leading order interactions and mass effects have 



been evaluated in (Philipsen and Zeidlewicz, 20101 



2.5 The quenched limit of QCD and Z(7V)-symmetry 

In the quenched limit, i.e. when quarks arc infinitely heavy, — > oo, the QCD 
partition function reduces to that of a pure gauge theory plus static quark fields. Let 
us examine the consequences of the compact temporal direction for gauge symmetry. 
The action is of course invariant under standard gauge transformations, 

S,[U^]^Sg[U] with U^^ix) = g{x)U^{x)g-\x + fi), g{x) e SU{N), (2.29) 

with our earlier periodic boundary conditions 

[/^(T,x)==[/^(T + iV„x), .g(r,x) =.g(r + iV„x) . (2.30) 

However, with the temporal boundary in place, we can also consider gauge transforma- 
tions with topologically non-trivial matrices g'{x), i.e. matrices that cannot be taken 
to unity by a smooth change of the parameters characterising the group elements. 
Consider a boundary condition on the transformation matrices, which is only periodic 
up to a constant matrix h, 

g'{T + Nr,x) = hg'{T,x) , h € SU{N) . (2.31) 

Such a g'{x) does not go into itself when winding around the torus once, but picks 
up a "twist" factor h. After a gauge transformation with g', the gauge links behave 
across the boundary as 

Uli'{T + Nr,x.) = hUf^{N^,^)h-^ . (2.32) 

They are only consistent with the periodicity requirement if \h, Uf ] = 0. This is 
satisfied if h is proportional to the unit matrix, i.e. it is in the centre of the group, 

27r7z 

h = zleZ{N), z^cxpi — , n € {0,1,2,. ..A^- 1} . (2.33) 

Thus, pure gauge theory at finite temperature is invariant under gauge transformations 
with non-trivial winding through the temporal boundary, for any global twist factor in 
Z{N). Note that this is not a symmetry of the QCD Hamiltonian. The latter is acting 
on a particular time slice and does not know about temporal boundary conditions. The 
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Hamiltonian generates translations in Euclidean time and its symmetries are defined 
on the space orthogonal to that. Rather, the centre symmetry discussed here is a 
symmetry of the space- wise Hamiltonian H^, Eq. (|2.8p . which acts on a z-slicc and 
thus is sensitive to the boundary in the temporal direction. 

Since the centre symmetry is related to non-trivial gauge transformations winding 
through the temporal boundary, gauge invariant observables can only be sensitive to it 
if they wind too. Such an observable is a Wilson line in the temporal direction closing 
onto itself, a Polyakov loop, 

Liic)^Y[Uo{x) . (2.34) 

Xq 

Physically, it corresponds to the propagator of a static quark. Under gauge transfor- 
mations, 

L^^) = gix)Li^)g-\x), 

Lf'(x) g'(l,x)i(x)g'"\l-f iV,,x) =5'(l,x)L(x).g'"'(l,x)/i-i . (2.35) 

Thus we see that the traced Polyakov loop is gauge invariant under topologically trivial 
gauge transformations, while it picks up a centre element when transformed with a 
winding transformation, 

TtL^ = TtL, TtL^' = z*TtL . (2.36) 

The Polyakov loop emerges naturally in the QCD path integral to leading order in the 
hopping expansion. For example, the partition function for a pure gauge theory with 
a static quark sitting at x is 

Zq= J DU Tri(x) e-^«[^l . (2.37) 



Hence, 



(Tri) = - I DU TrL e-^« = ^ = c-iFQ-Fo)/T ^ 
Z J Z 



(2.38) 



and the expectation value of the Polyakov loop gives the free energy difference between 
a Yang-Mills plasma with and without the static quark (McLerran and Svetitsky, 
1981). From this we can immediately infer two limiting cases: For T — > Yang-Mills 
theory is confining and it would cost infinite energy to remove the quark to infinity, 
i.e. Fq = oo and therefore (TrL) = 0. On the other hand, T — > oo corresponds to 
/3 — >■ oo, for which Uq ^ \ and (TrL) — >■ Trl = N . Clearly, a non-zero expectation 
value is no longer invariant under centre transformations and signals the spontaneous 
breaking of centre symmetry. Therefore, QCD in the quenched limit has a true (non- 
analytic) deconfinement phase transition corresponding to the breaking of the global 
centre symmetry, and the average of the Polyakov loop is the corresponding order 
parameter. 

If we add dynamical quark fields to the theory, they will behave as 

V'»(x) = g(x)V'(x), V(T + A^r,x) = -V(t-,x), V^'(T + ^r,x) -/#(r,x) . (2.39) 

Since statstical mechanics requires anti-periodic boundary conditions for fcrmions, the 
trivial h — 1 is the only permissible choice, i.e. there is no centre symmetry in the 
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presence of dynamical quarks. Physically, if there are dynamical quarks, their pair 
production screens the confining force (it leads to string breaking) and Fq is finite. 
Correspondingly, (TrL) 7^ for all temperatures and the Polyakov loop is no longer 
a true order parameter. In this case, a non-analytic phase transition as a function of 
temperature is not necessary, confined and deconfined regions may also be analytically 
connected by a smooth crossover. 



2.6 The chiral limit 

Chiral symmetry and its breaking on the lattice has been the subject of another 
set of lectures at this school. Let us merely summarise what we need for the finite 
temperature discussion. In the limit of zero quark masses the classical QCD Lagrangian 
in the continuum is invariant under global chiral symmetry transformations, the total 
symmetry being [/^(l) xS'/7l(A^/) xS'C/ij(-/V/). The axial J7^(l) is anomalous, quantum 
corrections break it down to Z{Nf). The remainder gets spontaneously broken to the 
diagonal subgroup, SUL{Nf) x SUB.{Nf) — >• SUv{Nf), giving rise to Nj — 1 massless 
Goldstone bosons, the pious. The order parameter signalling chiral symmetry is the 
chiral condensate. 

It is nonzero for T < Tc, when chiral symmetry is spontaneously broken, and zero 
for T > Tc- Hence there is a non- analytic finite temperature phase transition corre- 
sponding to chiral symmetry restoration. For non-zero quark masses, chiral symmetry 
is broken explicitly and the chiral condensate {4'tp) 7^ for all temperatures. Again, 
in this case there is no need for a non-analytic phase transition. Of course, for most 
lattice fermions chiral symmetry is either reduced (staggered) or broken completely 
(Wilson), rendering the behaviour in the chiral limit a most difficult subject of study. 



2.7 Physical QCD 

QCD with physical quark masses obviously does not correspond to either the chiral 
or quenched limit. The Z{3) symmetry as well as the chiral symmetry are explicitly 
broken. Nevertheless, physical QCD displays confinement as well as three very light 
pions as "remnants" of those symmetries. In the presence of mass terms there is no 
true order parameter, i.e. the expectation values of the Polyakov loop as well as the 
chiral condensate are non-zero everywhere. Hence the deconfined or chirally symmetric 
phase is analytically connected with the confined or chirally broken phase, and there is 
no need for a non-analytic phase transition. The following questions then arise, which 
should be answered by numerical simulations: for which parameter values of QCD is 
there a true phase transition, and what is its order? Are confinement and chirality 
changing across the same single transition or are there different transitions? If there is 
just one transition, which is the driving mechanism? If there is only a smooth crossover, 
how do the properties of matter change in the different regions? 
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2.8 Strong coupling expansions at finite T 

Strong coupling expansions are well known from spin models and QCD at zero temper- 
ature. In contrast to the asymptotic series obtained by weak coupling expansions, they 
yield convergent series in the lattice gauge coupling /3 = 2N/g'^ within a finite radius 
of convergence. The series approximate the true answer the better the lower (3, i.e. for 
a fixed Nr the lower the temperature. In pure gauge theory, the convergence radius 
is bounded from above by the critical coupling of the deconfinement transition, /3c. 
Hence, strong coupling series are analytical low temperature results, complementary 
to weak coupling perturbation theory which is valid at high temperatures. 

A detailed introduction to strong coupling methods in the vacuum can be found 
in ( [Montvay and Miinster, 1994[ ). Here we merely summarise the main formulae for 
Yang-Mills theory and discuss the modifications for finite temperature applications 
( [Langelage, Miinster and Philipsen, 2008[ ). 

The Wilson gauge action can be written as a sum over all single plaquette actions, 

SM^Y. E /3(l-^RcTrC/p] ^^5,. (2.41) 

Note that analytically we are able to consider an infinite spatial volume, A^^ — > oo. 
The plaquette action has an expansion in terms of characters Xr(^^) = TrL'r(?/) of the 
representation matrices Dr{U) of the group elements U, 

exp-5p = co(/3)[l + ;^d.c,(/3)x.(C/p)] , (2.42) 

r^O 

and dr is the dimension of the representation r. With these ingredients the free energy 
density can be written as 

/^-^lnZ=-61nco(/3)-i ^ a(C) [] <i>(^.)"" • (2-43) 

C=(X"') i 

where il. = V ■ N^. is the lattice volume and cg is the expansion coefficient of the 
trivial representation, which has been factored out. The combinatorial factor a{C) 
is introduced via a moment-cumulant-formalism, and equals 1 for clusters C which 
consist of only one graph or so-called polymer Xi . The contribution of a graph Xi is 




(2.44) 



The quantity in equation (j2.43p is customarily called a free energy, even at zero phys- 
ical temperature, because the path integral corresponds to a partition function if one 
formally identifies the lattice coupling /? with 1/T. One may use the fundamental rep- 
resentations of the gauge groups, c/ = u, as the effective expansion parameter, which 
together with some higher ones can be expressed as series in the lattice coupling 

SU(2): u=^ + 0(/32) SU(3): u^^ + 0{l3^). (2.45) 
4 18 

Here we are interested in a physical temperature T = 1/ (gNt), realised by compact- 
ifying the temporal extension of the lattice. The physical free energy is then obtained 
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Fig. 2.1 Graph Xi contributing to tlie lowest order f{NT,u) of tlie expansion of tiie piiysical 
free energy density at finite temperature. 

by subtracting the formal (Nr ~ oo) free energy, which removes the divergent vacuum 
energy as in the continuum. Thus the physical free energy density reads 

f{Nr, u) = f{N^, u) - /(oo, u) . (2.46) 

The coirtributing polymers Xi have to be objects with a closed surface, since 

dU Xr{U) = 5r,0 . (2.47) 



This means the group integration projects out the trivial representation at each link. 
To calculate the group integrals one uses the integration formula 

dU XriUV)xriWU-') = ^XriVW) . (2.48) 

dr 

Because of the difference in Eq. (|2.46p , those graphs contributing in the same way to 
f{Nr) and / (oo) drop out of the physical free energy. This is true for all polymers with 
time extent less than N-r- The calculation thus reduces to graphs with a temporal size 
of Nr on the finite N-r lattice, and graphs spanning or extending N^ on the infinite 
lattice. Such graphs contribute either to f{Nr) or to /(oo) (and in some cases to both). 
It is therefore clear that the strong coupling series for the physical free energy starts 
at a higher order than the formal zero temperature free energy. Moreover, the order 
of the leading contribution depends on Nt- 

The lowest order graph existing due to the boundary condition on the finite A',- 
lattice, but not on the infinite lattice, is a tube of length Nt with a cross-section of 
one single plaquette, as shown in Fig. l2.1l It forms a closed torus through the periodic 
boundary and thus gives a non-vanishing contribution, which is easily calculated to 
be $(Xi) = u^^^ . We need to sum up all such graphs on the lattice. There are three 
spatial directions for the cross section of the tube, giving a factor of 3. Translations 
in time take the graph into itself and do not give a new contribution, while we get 
V^{Xi) from all spatial translations. Together with the 1/fi in Eq. (j2.43p this gives 
a factor of 1/Nt- The contribution of all tubes with all plaquettes in the fundamental 
representation for SU (2) is thus 

^{X,)^^u^^^, (2.49) 
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which is - up to a sign - also the leading order result for the physical free energy. For 
SU (N) with > 3 we have an additional factor of 2 because there are also complex 
conjugate fundamental representations. Thus, in the strong coupling limit m — > the 
free energy density and pressure are zero. 

The leading correction comes from tubes with inner plaqucttes, higher orders have 
local decorations of additional plaqucttes either in the fundamental or in higher repre- 
sentations. For the interesting case of SU{3), these contributions up to the calculated 



orders are (Langelage, Miinster and Philipsen, 2007) 



Nr 
Nr 



2048 2048 
256 



, (2.50) 



with 6 = 1 - 3u - 6z; + 8w, c = 1 + 3u + 6w + 8w - 18^^ and v = (5^/452 + 0{l3'^),w = 
/3^/288 + 0{(3^). This series is valid only for Nt > 5, for smaller Nt there are modifi- 
cations coming from polymers with cross-sections larger than one plaquette. 

2.9 The strong coupling regime as an ideal hadron gas 

It is now interesting to ask how the QCD pressure can be interpreted in the strong 
coupling regime. From the Wilson action it is clear that the strong coupling limit is 
also non-interacting. However, as we have noted already, in this limit the pressure is 
zero. Considering strong but finite couplings, let us recall the first orders of the T = 
glucball mass calculations ( [Miinster, 1981[ |Seo, 198^ , 

o ^ 2 27 3 ^ 4 297 5 858827 g 47641149 ^ 

m[E--) ^ -41n u - 3. + 9.^ - ^ _ ]^^^. ^ 1104587 29577789 

^ ' 2 2 10240 71680 

/^+-^ o 9 3 98 4 33 5 36771 g 117897 7 ,^ 

m(T,+ ) = -4 In u + 3?^ + -u'^ iC ^ ^ u' , (2.51) 

^ ^ ' 2 4 4 1280 448 ' ^ ' 

where the arguments on the left side denote the representations of the point group, 
which map into the different spin states ( [Montvay and Miinster, 1994[ ). We observe 
that the expansion of the free energy can be written 

f{Nr,u) = - — \e-"^{At + )N^ ^^^-m^E+^)N^ ^^^-ra{Tr)NA {l + 0{u^)) . (2.52) 

The prefactors before the exponentials correspond to the number of polarisations of the 
respective glueball states. Note that higher spin states start with ~ 6 In u ( [Schor, 1983[ ), 
thus contributing to the order ^ u^^* or higher in the free energy. Hence, through two 
non-trivial orders our result is that of a free glucball gas, modified by higher order 
corrections. By employing a leading order hopping expansion, the same conclusion can 



be reached when heavy quarks are added (Langelage and Philipsen, 201061. This is 
a rather remarkable result. It allows to sec from a first principle calculation that the 
pressure is exponentially small in the confined phase, and that it is well approximated 
by an ideal gas of quasi-particles which correspond to the T ~ hadron excitations. 



3 

Some applications 



3.1 The equation of state 

Energy density e(r) and pressure p{T) as a function of temperature are certainly 
among the most fundamental thermodynamic quantities of QCD governing the ex- 
pansion of the plasma in the early universe as well as in heavy ion collisions. Let us 
use the ideal gas results to develop some intuition about what will happen in QCD. 
In the high temperature limit T — >■ oo, we have a gas of no n- interacting gluons and 
quarks with 

P ^ (lQ+ll2Nt^ — . (3.1) 

On the other hand, as T — >■ we consider a hadron resonance gas model. For temper- 
atures significantly below ^ 200 MeV, only the pions are relativistic, which come in 
three charge states with spin zero, v = 3. A gas of non-interacting pions has pressure 

4 = 3-. (3.2) 

Thus, for QCD we expect the pressure to change as a function of temperature from a 
small to a large value as a signal of deconfinement. 

For the fully interacting case, we need to compute the free energy density and 
the pressure from the QCD partition function. A technical obstacle here is that, in a 
Monte Carlo simulation, one cannot compute the partition function directly, since all 
expectation values are normalised to Z, 

(O) = Z-^Tj:{pO) . (3.3) 

So we need to specify some observable O to access the partition function. The most 
frequently used d etour is called the integral method ( Engels, Fingberg, Karsch, Miller 
and Weber, 1990), in which a derivative of the free energy is calculated and then 
integrated, 

4' =-ird.^^!M(M. (3.4) 



On the lattice, it is convenient to take the derivative with respect to the gauge coupling 
instead of temperature. 
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Fig. 3.1 Flavour dependence of the pressure for A^^ ~ 4 lattices compared to a continuum 
extrapolated pure gauge result. From ( |Karscti, Laermann and Peikert, 2000[ ). 

= -N^ d/3' (3(Tr[/^ + TrC/;) -6(Tr[/p)T=o) . (3.5) 

Now wc simply need to measure expectation values of temporal (Up) and spatial (Up) 
plaqucttcs with sufficient accuracy, so they can be integrated numerically. Note that 
this introduces a lower integration constant, which needs to be fixed for the result 
to be meaningful. While we do not know /(/3o) from first principles, we can choose 
/3o corresponding to a temperature below the phase transition, where the free energy 
should be well modelled by a weakly interacting hadron gas. For r<130 MeV, all 
hadrons become non-relativistic and we can approximate the pressure by zero. Note 
however, that this procedure is not good for the chiral limit, when pions become 
massless and stay relativistic down to very low temperatures. 

Another difficulty is that strong discretisation effects are to be expected. At high 
temperature the relevant partonic degrees of freedom have momenta of order ttT ~ 
TT/{aNt) on the scale of the lattice spacing, by which they are strongly affected. For 
the equation of state it is therefore particularly important to gain control over these 
effects and carry out the continuum limit a — > 0. This motivates the use of improved 
actions, designed to minimise cut-off effects in the approach to the continuum. 

The results of a computation of the pressure with an improved action (Karsch, 
Laermann and Peikert, 2000) are shown in F'ig. 13.11 The data have been obtained for 
Nf = 2,3 with (bare) mass niq/T = 0.4 as well as for Nf = 2 + 1 with a heavier 
mass rUg/T = 1 on a very coarse lattice, Nt = 4. Nevertheless, interesting qualitative 
features can be observed. For comparison, continuum extrapolated pure gauge results 
are also included. The figure shows a rapid rise of the pressure in a narrow transition 
region. The critical temperature as well as the magnitude of p/T* refiect the number 
of degrees of freedom liberated at the transition. This last conclusion is firm, since 
the pressure also rises for fixed temperature when light quarks arc added to the the- 
ory, consistent with the behaviour in the Stefan-Boltzmann limit. Another interesting 
feature is that the curves fall short of the ideal gas values, i.e. interactions are still 
strong just above Tc- An important question then is whether these features survive in 
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the continuum limit. In pure gauge theory this can be firmly established by numerical 
extrapolation while for light dynamical quarks this is the subject of ongoing simulation 
programs. 

3.2 Screening masses 

Essentially all static equilibrium properties of a thermal quantum field theory are en- 
coded in its equal time correlation functions. These are quantities that are well defined 
and calculable to good precision by lattice methods. Unfortunately, these quantities 
are not directly accessible in heavy ion collision experiments. Nevertheless, their theo- 
retical knowledge provides us with the relevant length scales in the plasma, from which 
conclusions about the active degrees of freedom and their dynamics may be drawn. 

The concept of screening is most easily introduced in a QED plasma. Suppose we 
insert a static external charge in a plasma of freely moving charges. These will be 
attracted or repelled, arranging themselves such that the polarisation cancels out the 
field of the static charge, which therefore is screened. Beyond a certain distance called 
a screening length, a particle in the plasma will not feel the presence of the external 
charge. The screening length is defined by the spatial exponential decay of the equal 
time correlator of the electric field, 

lim {Ei{x,t)Ej{0,t)) = const, e"™"!^!, (3.6) 

x^oo 

and the inverse screening length m^i is the Debye mass. In perturbation theory, niD 
corresponds to the pole of the propagator at Matsubara frequency uj„ = 0, i.e. it 
corresponds to the electric gluon mass, mjj = niE, Eq. (|1.15|) . Its Fourier transforma- 
tion leads to the Debye-screened potential of a static charge, 

^(^^ = (^p^+n„o(o,p) = ^ • (^-'^ 

The same concept has been carried over to QCD where it is speculated to be re- 
sponsible for the loss of confinement for heavy quarkonium states by screening of colour 
charges in the plasma ( Matsui and Satz, 1986[ ). However, there are some conceptual 



difficulties with the definition of the Debye mass in QCD. First, colour-electric fields 
are no physical observables since they are gauge- dependent, and so is their correlator. 
One might argue that the pole mass of this correlator is still gauge-invariant order by 
order in perturbation theory and consider the perturbative series ( |Rebhan, 1993D , 

mn = rn^n"" + f ff^^ln ^ + cg^T + 0{g') . (3.8) 
47r 

In this case one encounters the Linde problem already in next-to-leading order: the 
coefficiet c receives contributions from all loop orders and cannot be evaluated in 
perturbation theory. 

For a non-perturbative evaluation of screening, one therefore generalises the con- 
cept to gauge-invariant sources, averaged over Euclidean time, 

<5W = ]^ E (O(x)O(O)), = ^ c„ e-^^-'l-l . (3.9) 
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Fig. 3.2 Left: Leading grapli contributing to tiie plaquette correlation function. Right: Graph 
contributing to tlie finite temperature effect of the correlator. 

These are correlations in space, and the connected parts therefore fall off with the 
eigenvalues of the space- wise Hamiltonian H^, Eq. (|2.8I) . In analogy to QED, these 
masses correspond to the inverse length scale over which the equilibrated medium is 
sensitive to the insertion of a static source carrying the quantum numbers of O. Because 
of the compact Euclidean time direction at T > 0, the continuum rotation symmetry 
of the hypertorus orthogonal to the correlation direction is broken down from 0(3) 
to 0(2) X Z{2), corresponding to rotations in the (x, j/)-plane and reflections of r. 
Correspondingly, screening masses are classified by quantum numbers J^'" , where 
J, P, C are standard spin, parity and charge conjugation, while R is associated with 
the Euclidean time reflection. The appropriate subgroup for the lattice theory is x 
Z{2) = Df^. The irreducible representations and the classification of operators have 



been worked out for pure gauge theory (Grossman, Gupta, Heller and Karsch, 1994 



Datta and Gupta, 1998[ ) as well as for staggered quarks ( [Gupta, 1999 1 . 

Just as in T = spectrum calculations, we can consider either glueball or me- 
son and baryon like operators. Once again it is useful to develop some intuition by 
considering the limits of low and asymptotically high temperature. In the latter case, 
standard perturbation theory should apply and we can evaluate the correlators using 
quark and gluon propagators. For a mesonic operator the exponential decay at large 
distances is then dominated by two quark lines, with an effective mass given by their 
lowest Matsubara frequencies, ~ ttT. Neglecting bare quark masses, a meson correla- 
tor in the non-interacting limit then decays with Mqq 2'kT. Perturbative corrections 
lift the degeneracy and appear to generally shift this towards larger values (Laine and 
Vepsalainen, 2004). For purely gluonic operators, the lowest Matsubara frequency is 
zero and the exponential decay is determined by dynamically generated mass scales. 
For operators like TrAg the propagators feature the perturbative Debye mass, and 
hence Mj^2 2mE- However, for operators involving Ai, we arc faced with the Linde 
problem and a perturbative evaluation is not possible. 

On the other hand, at low temperatures we have small /3 and a calculation by 
strong coupling methods is feasible. As an example, let us consider the colour-electric 
field correlator ( TrPg'', (x) TrFg" (y) ) , which is in the Jf*^ = Oj+ channel (T denotes 
reflection in Euclidean time) containing the ground state and the mass gap. On the 
lattice, this corresponds to a correlation of temporal plaquettes, and the quantum 
numbers under the point group Df^ are Af^. Temporarily assigning separate gauge 
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couplings to all plaquettes, the correlator can be defined as 



2 



OPiOp2 



(3.10) 



At zero temperature the exponential decay is the same as for correlations in the time 
direction, and thus determined by the giueball masses, the lowest of which may be 
extracted as 

m = - lim - \n C{z). (3.11) 

z-S-oo z 

The leading order graphs for the strong coupling scries at zero temperature arc shown 
in Fig. 13.21 (left'). This leads to the lowest order contribution: 

C(z) = Au'^^ = .4e-"=^ (3.12) 

Thus, to leading order for the giueball mass is nis — —4 In u{/3). 

Now we switch on a physical temperature, i.e. keep the lattice volume compact in 
the time direction. As in the case of the free energy, we are here only interested in the 
temperature effects, i.e. in the mass difference 

Am{T) = Tn{T) - m{0) = - lim -[\nC{T; z) - lnC{0; z)] (3.13) 

z^oo z 



= — lim 

2— >00 



1,M , AC(T;.) 



C(0;z) 



(3.14) 



with AC(T;z) = C{T;z) — C(0;z). A typical graph contributing in lowest order to 
this difference is shown in Fig. 13.21 (right). Summing up all leading order graphs gives 
( Langelage, Miinster and Philipsen, 2008[ ) 



Am(r) = (3.15) 

i.e. the screening masses decrease compared to their T = values. Again, the leading 
order result is generic for all SU{N) and quantum number channels. As in the case of 
the free energy, the difference only receives A^T--dependent higher orders of u, leading 
to weak temperature effects. We conclude that in the confinement phase the lowest 
screening masses in each quantum number channel should be close to the corresponding 
zero temperature particle masses, with a significant temperature dependence showing 
up only near Tc- This explains the findings of numerical investigations of the lowest 
screening mass in SU{3) gauge theory, which for temperatures as high as T = 0.97Tc 
see very little temperature dependence, Am(T)/m(0) > 0.83 ( |Datta and Gupta, 1998D . 

Away from the high and the low temperature limits, numerical simulations have to 
be performed. One particularly interesting aspect of screening masses is that they per- 
mit to study chiral symmetry restoration across the quark hadron transition by looking 
at degeneracy patterns. For example, consider the lowest lying scalar, pseudo-scalar, 
vector and axial-vector mesons. Chiral SU{2) x SU{2) implies degeneracy between the 
scalar and vector, while an intact f/A(l) implies degeneracy between the parity flipped 
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states. A calculation of these masses for temperatures around Tc will thus allow insight 
into the details of chiral symmetry restoration. 
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3.3 The free energy of a static quark anti-quark pair 

Similar to a single static quark, Eq. (|2.37p and the discussion following it, we can 
consider a static quark anti-quark pair in a Yang-Mills plasma, 

(TrLT(x)Tri(O)) = = cxp ^ . (3.16) 

(In the following we drop the subtraction of Fq from the notation, which is always 
implied). Such a system is particularly interesting because it represents the non- 
relativistic limit of heavy quarkonia in the plasma, which are routinely produced in 
heavy ion collisions. In order to clarify the meaning of this quantity, let us consider 
its spectral decomposition. We use again the transfer matrix formalism, but with a 
slight modification accounting for the static sources. Let us fix to temporal gauge, 
C/o(t, x) = 1,T = 1, ...A^^ — 1, which we are allowed to do on all temporal links but 
those in one time slice. On the gauge fixed time slices, we have the Kogut-Susskind 
Hamiltonan Hq, which acts on the Hilbert space of states with static sources, from 
Eq. 

(To),+l,, EEC-"^^" =exp-i[[/,;(T+l),l,[/,(T)] . (3.17) 
Now we can be rewrite the Polyakov loop correlator exactly as 

(TrLl'(x)TrL(O)) = (3.18) 

iTr(^ro^-^| DUo{N^) C/oL(^.,x)C/o,,(^r,0)e-^[^^(i)-^«(^^)^^'(^^)l) 

Next, we employ gauge invariance of the kernel of the transfer matrix, Eq. (j2.3|) . under 
gauge transformation in the upper timeslice, 

Rig)L = i[[/f (r + 1), C/o(T)g-\ [/.(r)] = L[[/,(t + 1), C/o(t), [/,(t)], (3.19) 

where R{g)\ip) = imposes a gauge tranformation on the states it acts on. Choosing 
g(x) = Uo{Nt,x.), we then obtain 

(TrLt(x)TrL(O)) = |tr(ro^-P„„^/3) , (3.20) 



where we have introduced the projection operator 
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Po.p^.. = J Dg gl^i^)g^,iO)Rig) . (3.21) 

This operator annihilates all wave functions not transforming as 

i;0^[U^]^g{^)p,^P,s[U]gliO). (3.22) 

Specifically, 

Pa^fiuHys) = J^Sfjjdfj^sliJa,^) ■ (3.23) 

Thus P projects onto the subspace of states with a colour triplet sitting at x and an 
anti-triplet at 0. Inserting a complete set of eigenstates of the Hamiltonian Hq, we 
then find 

(TrLt(x)TrL(O)) = ^ E Mn^c,) e"^^ = ^E^~^^ ' (^.24) 

We easily recognise this as a ratio of partition functions, the numerator being just the 
Boltzmann sum over all energy levels £^^'^ , which are eigenvalues of Hq in the presence 
of the QQ pair. Let us now consider the zero temperature limit of this expression. 



(3.25) 

Here we have identified the static potential with the ground state energy of a static 
quark anti-quark pair above the vacuum, F(|x|) = £^^'^(|x|) — Eq. This is precisely 
the same quantity which at zero temperature determines the ground state of the 
exponential fall-off of the Wilson loop. Thus, we can extract the same quantity from 
a Polyakov loop correlator, Eq. (|3.24p . in the limit Nr — oo. At finite temperatures 
instead, the Polyakov loop correlator corresponds to the Boltzmann- weighted sum over 
all excited states of the static potential, and hence to a free energy. This free energy 
is often called a T-dependent potential, V{r,T) = FQQ{r,T), even though this is not 
quite justified, as we will discuss below. 

The Polyakov loop correlator is readily simulated, with results as in Fig. 13.31 It 
gives a linearly rising free energy in the confined phase, whose effective string tension 
reduces with temperature, while in the deconfined phase the potential is screened, thus 
exhibiting clearly the phase transition from a confined to deconfined phase. However, 
the screened free energy does not correspond to the Debye-screened potential from 
Eq. p.7p . as becomes apparent when considering its spatial decay at high T. Fitting 
the screened free energy to 

EQQ^_^^-miT)r (3 26) 

gives d 1.5 and m = i.e. the screening mass corresponds to the lightest 

glueball channel ( [Hart, Laine and Philipsen, 2000[ ). This can already be seen in per- 
turbation theory, where the leading term is by two-gluon exchange and thus m = 2tod. 
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Fig. 3.3 Static quark anti-quarlc free energy/potential, Eq. (|3.16|l . for T < Tc (left) and 
T > Tc (riglit), on 32 X 4 witli Symanzilt improved action (Kaczmarek, Karscii, Laermann 
and Liitgemeier, 2000). 

Non-perturbatively it follows from the fact that the traced Polyakov loop is a gauge- 
invariant operator coupling to all powers of gauge fields, hence its correlator at large 
distance is dominated by the lightest gauge-invariant gluonic screening mass. 

For the following, let us consider explicitly = 3 colours. In order to find a 
thermal potential in analogy to QED Debye screening, attempts have been made to 
decompose the traced and gauge invariant Polyakov loop into colour singlet and octet 
configurations of the static sources, 3x3 = 1 + 8, ( [McLerran and Svetitsky, 1981[ 



Nadkarni, 1986) 



"9 ^9 



-Fi{r,T)/T _ 



-Fg{r,T)/T 



i(TrLt(x)L(y)), 



3 

l(TrLt(x)TrL(y)) - l(TrLt(x)L(y)). 



(3.27) 



(3.28) 



Note that the correlators in the singlet and octet channels are gauge dependent, and 
the colour decomposition only holds perturbatively in a fixed gauge, which has moti- 
vated many gauge fixed lattice simulations. However, both options are unphysical at 
a non-perturbative level. To understand this, let us start from something physical and 
consider a meson operator in an octet state, 0° = '0(x)[/(x, Xo)T°i7(xo, y)?/'(y), with 
So the meson's center of mass. In the plasma the colour charge can always be neu- 
tralised by a gluon. In the correlators for the singlet and octet operators, we integrate 
out the heavy quarks, replacing them by Wilson lines. 



(O(x,y;0)Ot(x,y;7V^)) cx (TrLt(x)C/(x, y; 0)L(y)C/t (x, y; TV,)), 

1 



(O"(x,y;0)O"t(x,y;7V.)> « 



(TrLt(x)Tri(y)) 



(3.29) 



-^(TrLt(x);7(x,y;0)L(y)C/t(x,y;iV,)) 
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We have now arrived at gauge invariant expressions, because we used a gauge string be- 
tween the sources. The singlet correlator corresponds to a periodic Wilson loop which 
wraps around the boundary. The connection to the gauge fixed correlators is readily cs- 
tabhshcd, replacing the gauge string by gauge fixing functions, J7(x, y) = g~^{x)g{y). 
Thus, in axial gauge, J7(x, y) = 1 (and only there), the gauge fixed correlators are 
identical to the gauge invariant ones. 

Repeating the spectral analysis the full correlators take the form ( Jahn and Philipsen, 
20045 



-F.(..T)/T ^ ^1 ^(,,^|[/^,(x,y)[/t^(x,y)|n,„)e-^?''W/^ 



e-^^^^'^y^ = li^(„,,|C/^«,(x,y)C/I;(x,y)|n,„)e-<'''M/^ . (3.30) 

n 

The energy levels in the exponents are identically the same in Eqs. p.24l3.30p and 
correspond to the familiar gauge invariant static potential at zero temperature and 
its excitations. However, while Eq. p.24p is purely a sum of exponentials and thus 
a true free energy, the singlet and octet correlators contain matrix elements which 
do depend on the operators used, thus giving a path/gauge dependent weight to the 
exponentials contributing to Fi, Fg. Since the spectral information contained in the 
average and gauge fixed singlet and octet channels is the same, we must conclude 
that any difference between those correlators is entirely gauge dependent and thus 
unpliysical. 

The definition of a temperature-dependent static potential for bound states with 
the appropriate perturbative Debye screened limit has recently been achieved. How- 
ever, it requires a Wilson loop correlator in real time evaluated in a thermal ensemble 
( |Laine, Philipsen, Romatschke and Tassler, 2007 Brambilla, Ghiglieri, Vairo and Pe- 



treczky, 2008j~ Since this cannot be done by straightforward lattice simulations, we 
shall not further discuss it here. 

3.4 Phase transitions and phase diagrams 

A question of utmost importance for experimental programs ranging from heavy ion 
collisions to astro-particle physics are the different forms of nuclear matter for vari- 
ous conditions specified by temperature and baryon chemical potential, and whether 
these are separated by phase transitions. In statistical mechanics, phase transitions 
are defined as singularities, or non-analyticities, in the free energy as a function of its 
thermodynamic parameters. However, on finite volumes, free energies are always ana- 
lytic fun£tion£_aiid_tl^^ Lee and Yang, 



1952) states that singularities only develop in the thermodynamic limit of infinitely 
many particles, or V ^ 00. This is particularly obvious in the case of lattice QCD, 
whose partition function is a functional integral over a compact group with a bounded 
exponential as an integrand and without zero mass excitations. It is thus a perfectly 
analytic function of T, fi, V for any finite V. Hence, a theoretical establishment of a 
true phase transition requires finite size scaling (FSS) studies on a series of increasing 
and sufficiently large volumes to extrapolate to the thermodynamic limit. 
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Fig. 3.4 Conjectured phase diagram for Nf = 2 QCD witli finite ligfit quark masses. The 
physical case Nf = 2 + 1 is believed to qualitatively look the same, based on universality and 
continuity arguments as well as input from QCD-like models ( [Rajagopal and Wilczek, 2000[ |. 

Three different situations can emerge: a first order phase transition is characterised 
by coexistence of two phases, and hence a discontinuous jump of the order parameter 
(and other quantities), while a second order transition shows a continuous transition 
of the order parameter accompanied by a divergence of the correlation length and 
some other quantities, like the heat capacity. Finally, a marked change in the physical 
properties of a system may also occur without any non-analyticity of the free energy, in 
which case it is called an analytic crossover. A familiar system featuring all these possi- 
bilities is water, with a weakening first oder liquid-gas phase transition terminating in 
a critical endpoint with Z(2) universality, as well as a triple point where the first order 
liquid-gas and solid-liquid transitions meet. Similar structures are also conjectured to 
be present in the QCD phase diagram (Halasz, Jackson, Shrock, Stephanov and Ver- 



baarschot, 1998; Rajagopal and Wilczek, 2000[ ), Fig. 13.41 where asymptotic freedom 



suggests the existence of at least three different regimes: the usual hadronic matter, a 
quark gluon plasma at high T and low /i, as well as a colour super-conducting state 
at low T and high /i. 



3.5 The (pseudo-) critical coupling and temperature 

The first task when investigating a potential phase transition is to locate the phase 
boundary. Thus, for QCD with a fixed quark content, we are interested in the critical 
(or pseudo-critical) temperature where the transition from the confining regime to the 
plasma regime takes place. The method to locate a transition in statistical mechanics 
usually is to look for rapid changes of suitable observables 0(x) and peaks in their 
susceptibilities. Typical examples are the Polyakov loop, the chiral condensate or the 
plaquette, O S {TrL, 'ipij), TiUp, . . .}. Generalised susceptibilities for 0(x) in statistical 
mechanics are defined as the volume integral over the connected correlation function, 

Xo - 1 d'x ((O(x)O(O)) - (O(x))(O(0))) . (3.31) 

(Note that on the lattice this can be generalised to 4d. However, the thermal equilib- 
rium system likewise lives in three dimensions, which are the ones to be used for finite 
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Fig. 3.5 Polyakov loop (TrL) and chiral condensate {i)tp) together with their susceptibilities 
signal a transition in two flavour QCD. From ( jKarsch, 2002| . 

size scaling analyses in the following section) . If we instead consider the corresponding 
volume average, 

O = J d^x 0(x), (3.32) 

then because of its translation invariance the integration gives trivial volume factors 
and we obtain on the lattice 

Xo = - {Of) = NmSOf), (3.33) 

with the fluctuations SO = O — (O). At a phase transition fluctuations are maximal, 
hence the locations of the peaks of susceptibilities define (pseudo-) critical couplings, 
xiPc'Tif) = Xmax =^ l^cijnf), which can be turned into temperatures as discussed in 
Sec. 12.21 In practice, often the two-loop beta function is used as a short cut. although 
this becomes valid only when the lattice spacing is fine enough to be in the perturbative 
regime. 

Note that for an analytic crossover the pseudo-critical couplings defined from dif- 
ferent observables do not need to coincide. The partition function is analytic every- 
where and there is no uniquely specified "transition" . This holds in particular for 
pseudo-critical couplings extracted from finite lattices. As the thermodynamic limit is 
approached, the couplings defined in different ways will merge where there is a non- 
analytic phase transition, and stay separate in the case of a crossover, as illustrated 
in Fig. 13.61 for a putative QCD phase diagram. 

3.6 Universality, finite size scaling and signals for criticality 

A fascinating phenomenon in physics is the "universality" exhibited by physical sys- 
tems near critical points of second order phase transitions. It is due to the divergence 
of the correlation length, which implies that the entire system acts as a coherent 
collective. Hence, microscopic physics becomes unimportant, the collective behaviour 
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Fig. 3.6 Location of pseudo-critical parameters defined by different observables. In the 
infinite volume limit, the lines merge for a true phase transition, but stay separate for a 
crossover. 

is determined by the global symmetries and the number of dimensions of the system. 
With the divergence of the correlation length any characteristic length scale disappears 
from the problem, and thermodynamic observables in the critical region obey scale in- 
variant power laws. For example, at a second order ferromagnetic phase transition, 
the magnetization (the order parameter) vanishes as M ~ |t|^ with t = {T — Tc)/Tc, 
while the specific heat, the magnetic susceptibility and the correlation length diverge, 
C ~ X ~ 1^1"''' and ^ ^ respectively. The critical exponents a,/3,7, and 

similar ones for other quantities arc the same for all systems within a universality 
class. The latter are usually labeled by spin models, since those can be readily solved 
numerically. Unfortunately for most systems, and particularly QCD, it is not obvi- 
ous which global symmetries the system will exhibit at a critical point, since those 
are a dynamically determined subset of the total symmetry of the theory. Moreover, 
since there is no true order parameter in the case of finite quark masses, fields and 
parameters of QCD map into those of the effective model in a non-trivial way. 

Nevertheless, general scaling properties can be derived, we follow ( Karsch, Laer- 
mann and Schmidt, 2001). In the vicinity of a critical point, the dynamics of QCD 
will be governed by an effective Hamiltonian in analogy to a spin model, with energy- 
like and magnetisation-like (extensive) operators E, M which couple to two relevant 
couplings. 



T 



= tE + hM. 



(3.34) 



In the case of the Ising model r and h are proportional to temperature and mag- 
netic field, respectively. The power law behaviour for thermodynamic functions near 
a critical point follows from the scaling form of the singular part of the free energy. 



fs{T,h)^b-''Mb''-T,b''-h), 



(3.35) 



with a dimcnsionlcss scale factor b = LT = Ng/N^, the spatial dimension d and the 
dimensions of the scaling fields D-r, Dji. Using general scaling relations between those 
and the critical exponents. 



(3.36) 
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one derives the scaling of the susceptibilities as 



TdT^ 




(3.37) 



(Note the different volume factors between Eqs. p.33l3.37p . which is due to the fact 
that Eq. p.37p uses extensive variables while in Eq. p.33p we use averages.) In QCD, 
the operators E, M will be mixtures of the plaquette action and the chiral condensate 
(and possibly higher dimension operators), while r and h are functions of P^rrif^fj,. 
Thus, when measuring susceptibilities constructed from operators as in Eq. p.33p . 
these will in turn be mixtures of E and M. Therefore in the thermodynamic limit all 
of them will show identical ESS behaviour which is dominated by the larger oi a/v 
and "i/v- For the universality classes relevant for QCD, Z{2) and 0(4),0(2), this is 
l/v- 

By contrast, at a first order phase transition the fluctuations diverge proportional 
to the spatial volume. 



while for an analytic crossover the peaks of susceptibilities saturate at some finite 
maximal value in the thermodynamic limit. 

Another possibility to detect phase transitions is to evaluate the partition function 
for complex values of the coupling /3 and look for Lee- Yang zeroes of the partition 
function, Z{(i) = 0. On a finite volume, these will be at complex values of the coupling, 
since the partition function is analytic for real parameter values. On larger volumes a 
Lee- Yang zero approaches the real axis if there is a true phase transition. However, such 
evaluations require multi- histogram and reweighting techniques (cf. Sec l4.3p . where 
the values of the observable at different parameter values are calculated from one 
simulation point, a technique which has to be handled with care ( |Ejiri, 2d06| ). Fig. 13.71 
(left) shows an example of Lee- Yang zeroes for SU{3) pure gauge theory at complex 
values of the coupling. 

A more straightforward but numerically very expensive observable to determine 
the order of a phase transition is the Binder cumulant constructed from moments of 
fluctuations of the order parameter. 



It is to be evaluated on the phase boundary (3 = (^cirri f , fi) , where the third moment 
of the fluctuation vanishes, {{SO)^){(3c) = 0. This observable is particularly well suited 
to locate the change from a first order transition to a crossover regime as a func- 
tion of some parameter, like quark mass or chemical potential. In the infinite volume 
limit, ,64 —> 1 or 3 for a first order transition or crossover, respectively, whereas it 
approaches a value characteristic of the universality class at a critical point. For Z{2) 
(or 3d Ising) universality one has B4 — > 1.604. Hence, when examining the change of 
a phase transition from a weakening first order transition to a critical end point and 



(3.38) 



{{SOW ■ 



(3.39) 
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Fig. 3.7 Left: Lee- Yang zeroes in tlie complex /3-plane for SU(3) pure gauge tlieory (Ejiri, 
2006r 

Riglit: Tlie Binder cumulant as a function of some parameter x wiiicii takes tiie pliase 
transition from a first order to a crossover regime, passing through a critical point of 3d Ising 
universality. (E.g. x = m/ in the case of Nf — 3, fi = 0, cf. Fig. I3.8[l . 

a crossover as a function of x ~ anif^afi, is a non-analytic step function, which 
gets smoothed out to an analytic curve on finite volumes, with a slope increasing with 
volume to gradually approach the step function. Fig. 13.71 (right). Near a critical point 
the correlation length diverges as ^ ~ , where r is the distance to the critical point 
in the plane of temperature and magnetic field like variables. Since the gauge coupling 
is tuned to /3 = I3c{x), r ~ |a; — Xc\ and is a function of the dimensionless ratio L/^, 
or equivalently, {L/^Y/" . It can therefore be expanded as 

Bi{x,N,) ^hQ + bNll''{x^x'') + .... (3.40) 

Close to the thermodynamic limit all curves intersect at a critical i?4-value, and more- 
over V takes its universal value. Thus calculation of i?4 around a critical point gives 
access to two independent pieces of information governed by universality. 

3.7 The nature of the QCD phase transition for TV/ = 2 + i at ^ = 

The qualitative picture for the order of the QCD phase transition at zero baryon den- 
sity as a function of the quark masses is outlined schematically in Fig. 13.81 As discussed 
in Sees. 12.512.61 for Nf = 3 in the limits of zero and infinite quark masses (lower left 
and upper right corners), order parameters corresponding to the breaking of a symme- 
try can be defined, implying true phase transitions. One finds numerically at small and 
large quark masses that a first-order transition takes place at a finite temperature Tc- 
On the other hand, one observes an analytic crossover at intermediate quark masses. 
Hence, each corner must be surrounded by a region of first-order transition, bounded 
by a second-order line. The line bounding the chiral transitions is commonly called 
the chiral critical line, the one bounding the deconfinement transition is accordingly 
the deconfinement critical line. We know from simulations on Nr = 4 with staggered 
fermions ( |de Forcrand and Philipsen, 2007] ) that this line is to the lower left of the 
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Fig. 3.8 Order of the QCD phase transition as function of quark masses (m„_d,"is) at 
^I = 0. 



physical point, and simulations on Nr = 6 (de Forcrand, Kim and Philipsen, 2007 



[Endrodi, Fodor, Katz and Szabo, 2007[ ) show that it shrinks towards the lower left 
corner on finer lattices. A continuum extrapolation for the chiral critical line is not 
yet available. However, there are continuum extrapolated simulations with staggered 
quarks with physical masses that show that the /i = transition in QCD is a crossover 



(Aoki, Endrodi, Fodor, Katz and Szabo, 20061 



The situation is not yet settled for the chiral limit of the two flavour theory in the 
upper left corner. If the transition is second order, then chiral symmetry SU{2)l x 
SU{2)ii ~ 0(4) puts it in the universality class of 3d 0(4) spin models. In this case 
there must be a tri-critical strange quark mass m}"'^, where the second order chiral 
transition ends and the first order region begins. The exponents at such a tri-critical 



point would correspond to 3d mean field (Lawrie and Sarbach, 1984). On the other 



hand, a first order scenario for the chiral limit ofNf = 2 so far has not been conclusively 
ruled out. In fact, it has been shown in a model with the same chiral symmetries as 
Nf = 2 QCD and a tuneable C/4(l) anomaly that both scenarios are possible and 
the order of the transition depends on the strength of the anomaly at the critical 
temperature ( Chandrasekharan and Mehta, 2007] ) . 



4 

Lattice QCD at non-zero baryon 
density 



4.1 Implementing chemical potential 

So far wc have considered lattice QCD at finite temperature with net baryon number 
zero, i.e. with a complete balance between baryons and anti-baryons. This situation 
is reaUsed during the evolution of the early universe because the matter anti-matter 
asymmetry is so tiny there. However, there are interesting questions about the nature 
of dense nuclear matter in the centre of compact stars, and heavy ion collision exper- 
iments obviously also operate at finite baryon number. Thus, we now wish to use the 
grand canonical ensemble with a chemical potential /i for quark number Q (recall that 
each quark carries 1/3 of baryon number, Q = -B/3, /i = j-is/S), 



Tr e 



(fix ip(x)jo'>P{x) = I (fix 'ijj\x)'ip{x 



(4.1) 



For the following sections we recall that in Euclidean space time 7^ = 7;'^, {75, 7^i} = 0. 
Under charge conjugation, C = 7072, fermion number changes sign. 
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(4.2) 



Thus, /i > represents a net baryon number and /U < a net anti-baryon number. 
From this one derives an important symmetry of the partition function. Since the 
measure and the fermion action for /i = are invariant under charge conjugation, we 
find that the partition function is an even function of chemical potential. 



Z(/_i) = / DA^ Di^'-'DiP'-' exp - 



= J DA DipDip exp - 

= z{-^Ji) . 



l/T 



dxo Q 



Sg + Sf{p. = Q) + II 



1/T 



dxo Q 



(4.3) 



In a straightforward implementation of chemical potential on the lattice we would 
simply replace the integral by a lattice sum. 



(4.4) 



Implementing chemical potential 37 



However, when computing the energy density in lattice perturbation theory, one finds 
it to diverge in the continuum Umit, 

1 ^InZ^oo. (4.5) 



Vdi 



T) 



This happens after duly subtracting the divergent vacuum contribution, and thus the 
chemical potential term in Eq. (|4.4p appears to spoil renormalisability. This would vio- 
late our earlier observation that renormalisation of the vacuum theory is also sufficient 
for the theory at finite temperature and density. The reason is that the discretisation 
viola tes a generalised symmetry of the continuum theory ( Hasenfratz and Karsch, 
1983). Note that the term ^.Q in the continuum looks like the zero component of the 
fermion current, = ipj^ip, coupling to an external C/(l) gauge field, 

fiQ ~ —igAoio with Aq ~ i— . (4-6) 

9 

Chemical potential for quark number is equivalent to a classical electromagnetic field 
Aq with a constant imaginary value. The quark number term can therefore be absorbed 
into the covariant derivative of the Dirac action, which then is invariant under U{1) 
gauge transformations of the quark fields and Aq . This symmetry protects against new 
divergences and gets broken by a lattice implementation as in Eq. (|4.4|) . The problem 
is solved by implementing chemical potential in the same way as an external U{1) 
gauge field, namely as an additional temporal link variable 

Uo,e.t = c"3^« = c--^^ , (4.7) 

which multiplies all non-abelian temporal links. For example, the Wilson action with 
finite chemical potential then takes the form 

X 

-K [e''''V^(a;)(r - 7o)f/o(x)?/'(x - 6) + c-"'V(a; + 6)(r + 7o)t/J(a;)^(a;) 
■Y, [^(x)(r-7,)C/,(a;)^(a;+.7)+V^(a;+.7)(r + 7,)C/J(x)V(x)] j .(4.8) 



3 



As usual on the lattice, this discretisation is not unique, only the continuum limit is. 
For alternative ways of implementing chemical potential, see ( |Gavai, 1985[ ). It is now 
easy to verify that 

det(^([/t) +m + 7om) = det(^([/) + m - 70^) , (4.9) 

and since Sg[U'^] = Sg\U] and DC/t = DU we have Z{^) = Z{—i^i) on the lattice as 
well. 



38 Lattice QCD at non-zero baryon density 



4.2 The sign problem 

Let us consider what happens to the Dirac operator in the presence of chemical poten- 
tial. Reality of the fermion determinant follows from the 75-hermiticity of the Dirac 
operator, 

(^ + m)t =75(^ + m)75 . (4.10) 

This relation is satisfied in the continuum as well as for all standard lattice Dirac 
operators. Now consider a complex chemical potential, 

75(^ + m- 7oAi)75 = {-Ip + m + 7oA*) = {Ip + m + 70^^*)^ . (4.11) 

For the fermion determinant this implies, 

det(^ + m - 7oAi) = det*(^ + m + 7oAt*) , (4.12) 

which is complex unless Re/^ = 0. However, a complex fermion determinant cannot be 
interpreted as a probability measure for importance sampling and the evaluation of the 
path integral by Monte Carlo methods is spoiled. This is known as the "sign-problem" 
of QCD. 

It should be stressed that there is nothing wrong with the theory at finite /.i. In par- 
ticular, the partition function as well as the free energy and the other thermodynamic 
functions are all real after performing the path integral over the gauge fields, i.e. the 
imaginary parts of the fermion determinant in the background of gauge configurations 
average to zero. Being a property of the Dirac operator, the sign problem is generic 
for fermionic systems with particles and anti-particles or particles and holes, and is 
necessary for a correct description of the physics. For example, for the expectation 
value of the Polyakov loop we have 

(TrL) = e--^ = (ReTrL Redet Af - ImTrL Imdet Af)g , 
((TrL)*) = e--^ = (ReTrL Re det Af + ImTri Imdct A/)g , (4.13) 

where the angular brackets denote path integration with respect to the pure gauge 
action, (. ■ ■)g = J DU . . . exp — S'g[[/]. Thus, the free energy of a static quark or an 
anti-quark in a plasma differ if and only if Imdet 7^ 0, i.e. for finite chemical potential 
^. If there is no net quark or baryon number, it costs equal amounts of energy to 
insert a quark or an anti-quark into the plasma, whereas this is different if the plasma 
already has a net baryon number. The origin of the sign problem may thus be traced 
to the behaviour under charge conjugation, 

det(^ + m - 7om) det(^ + m + 70^^) . (4.14) 

The "sign problem" thus is only a problem for a Monte Carlo evaluation of the partition 
function with importance sampling methods. In some spin models, it can be solved by 
cluster algorithms, which are able to identify conjugate configurations so as to cancel 
imaginary parts during the simulation ( jWiese, 2002[ ). Unfortunately, these work only 
for special classes of Hamiltonians and QCD does not appear to be one of them. 



The sign problem 39 



Nf = 2 QCD at finite isospin density 



T 




quark-gluon 






plasma 




^ 9 
\ * 

\ 


9 




9 


9 — 




hadronic 


pion 




matter 


superfluid 



m„ Mi 

Fig. 4.1 Schematic phase diagram for QCD at finite isosipin density. 

Another idea is to simply avoid importance sampling and evaluate expectation values 
by means of stochastical quantisation, or Langevin algorithms (Karsch and Wyld, 
1985). However, also in this case there are difficulties, for complex actions there is no 
proof that the Langevin method converges to the correct expectation value. While the 
method seems to work for various model systems, there appears to be no successful 
application to QCD yet. 

Due to the symmetry properties of the determinant, there are a number of in- 
teresting special cases permitting a Monte Carlo treatment. As long as the fcrmion 
determinant is real, we arc able to simulate even numbers of degenerate flavours with 
it. Eq. ()4.12|) shows that this is the case for purely imaginary /i = i^i, fii S M, and we 
shall discuss later how this can be used to learn something about real ^. 

Another interesting situation for Nf = 2 is when the chemical potential for u and 
d quarks are opposite, i.e. /i„ = —fJ-d = M/j which corresponds to a chemical potential 
for isospin ( [Son and Stephanov, 200l| ). In this case we have for the determinants 



det(^ + TO - 7oM/) dct(^ + m + 70^/) = | det(^ + m - 7oM/)P > . (4.15) 



Note that isospin is not conserved in the real world due to the electroweak interactions. 
Nevertheless, one can investigate interesting finite density effects with this theory by 
actual simulations ( |Kogut and Sinclair, 2008"] ) . For example, one obtains a condensate 
of charged pions (tt") 7^ when > to^, Fig. 14.11 

Finally, the determinant is real for a two-colour version of QCD with gauge group 
SU{2). This is because this group has real representations, i.e. there is a matrix S 
such that ST'^S"^ = -T"* with S = cr^, the Pauh matrix. Now consider S = C-f^a'^, 
and using C^^C^^ — — 7J we have 
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S[lp + m - "f(3^]S ^ C75cr^[7^(9^ - igAp) + m - 7o/i](T^75C ^ 
= C'75[7^(9^ + igA*^) + m - 7om]75C"^ 

= [^ + m-7o/i1* . (4.16) 

For real ^ we have det M = det*M, likewise permitting simulations to study qualitative 
density effects ( [Hands, Sitch and SkuUerud, 2008[ ). 

In the following sections, the main concern will be with methods for the theory 
of immediate experimental interest, i.e. three-colour QCD with real quark number 
chemical potential. As will become apparent, all those methods side-step the sign 
problem by introducing approximations which work only for sufficiently small chemical 
potential, empirically one finds these methods to be good as long as < T. 

4.3 Reweighting 

Reweighting techniques are widely used in Monte Carlo simulations with importance 
sampling, in particular in order to interpolate data between simulation points when 



studying phase transitions (Ferrenberg and Swendsen, 19881. In the context of finite 



density physics, this technique is used to extrapolate to finite chemical potential. The 
basic idea is to rewrite the path integral exactly, 

Z{^i) ^ j DU dctA/(M) c-^'[^l ^ j DU detM(O) e"^"!^! 

^ ^\detM(0)/^^o ^ ' 

The angular brackets now denote averaging over an ensemble with the same measure 
DU det M(0) as for zero density, i.e. the path integral can be interpreted as an expec- 
tation value of a reweighting factor given by the ratio of determinants. While this is a 
mathematical identity, its practical evaluation by Monte Carlo methods turns out to 
be impressively difficult. First, a numerical evaluation of this path integral requires the 
calculation of the reweighting factor configuration by configuration and is expensive. 
This is aggravated by a signal to noise ratio which worsens exponentially with volume 
because the free energy is an extensive quantity F = V f, 

||^i.„,_£W_£<2).,,p_|,,M_,,„))^ (4,18) 

Hence the reweighting factor gets exponentially small as we increase V or ijl/T, re- 
quiring exponentially increased statistics for its determination. The problem can be 
illustrated with a one-dimensional Gaussian integral with a complex "action" , 

(4,19) 

Fig. 14.21 (left) shows the real part of the integrand for = and ^ = 30. It is obvious 
that the oscillatory bevahiour of the integrand will rapidly become prohibitive for 
Monte Carlo integration with growing /x. 
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Fig. 4.2 Left: The integrand of Eq. (|4.19p for different fi. Right: In order for reweighting to 
work, the probability distribution of the Monte Carlo ensemble must have sufficient "overlap" 
with the target integrand. 

The second problem is known as the overlap problem. When we evaluate the inte- 
gral with importance sampling, the integrand will be calculated the more frequently 
the larger its contribution to the integral. However, when we employ reweighting, the 
configurations will be generated with the measure DU detM(O), whose probability 
distribution is shifted compared to that of DU detM(/x), Fig. 14.21 (right). If we had 
infinite statistics, this would not matter as long as our algorithm is crgodic and the 
entire configuration space is covered. The point of Monte Carlo methods is, however, 
to get away with a few "important" configurations on which we evaluate our deter- 
minant and observables. But if the difference between the integrands at ^ = and 
H gets too large, the most frequent configurations are the unimportant ones. It is 
apparent that we only get a useful estimate for observables as long as there is sufficient 
overlap between the two integrands. The difficulty with this method is to realise when 
it fails. The statistical errors are determined by the fluctuations within the generated 
ensemble and will appear small, unless there are configurations that have tunneled far 
into the tail where the actual integrand is important. 

When searching for a phase transition at some fixed /3c(Ai), one-parameter reweight- 
ing uses an ensemble at /3c(/.t),/^ ~ 0, which is non-critical since Tc(/x) < Tc(0), thus 
missing important dynamics. Significant progress enabling finite density simulations 
was made a few years ago, by a generalisation to reweighting in two parameters (Fodor 
and Katz, 2002). The partition function is now rewritten as 



where the ensemble average is generated at ^ = and a lattice gauge coupling /3o, 
while a reweighting factor takes us to the values fj,, (3 of interest. In this way one can 
use an ensemble on the phase boundary at ^ = 0, which actually samples both phases, 
and then reweight along the (pseudo-) critical line to the desired value of /x. On physical 
grounds, this is drastically improving the overlap when attempting to describe a phase 
transition. 




(4.20) 
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Clearly, the choice of reweighting factors is not unique and one may ask for an 
optimal choice. A useful criterion is to minimise fluctuations in the reweighting factor, 



Eq. (|4.17p (|de Forcrand, Kim and Takaishi, 2003[). The optimal setup with a feasible 



Monte Carlo implementation is to split the determinant into modulus and phase, 
det M = I det M|e'^, and employ the modulus for the generation of the ensemble. Note 
that this is equivalent to an ensemble with finite isospin chemical potential. Thus, 

Z(^) = Z{^i,){c^'),, . (4.21) 

Numerical experiments show that for small to moderate chemical potentials, this choice 
is advantageous ( |de Forcrand, Kim and Takaishi, 2003[ ) compared to the standard pro- 
cedure. Moreover, it serves for interesting theoretical insights. There are indications in 
QC D (pplittorff, 2005[) and firm results in random matrix models ( Han and Stephanov, 
2008) that the phase transition to the pion condensate in the theory with finite isospin 
chemical potential is directly related to the strength of the sign problem, which be- 
comes maximal with an average sign near zero in the neighbourhood of that transition, 
thus signalling a breakdown of reweighting methods there. 

4.4 Finite density by Taylor expansion 

Another straightforward way to employ Monte Carlo ensembles generated at /i = to 
learn about finite ^ is by Taylor-expanding the partition function in ^/T . The grand 
canonical partition function is an analytic function of its parameters away from phase 
transitions. The Lee- Yang theorem tells us that on finite volumes there are no non- 
analytic phase transitions, and so the partition function must be analytic everywhere. 
It is then natural to consider power series m fi/T for thermodynamic obscrvables like 
the pressure ( |Allton et al, 2002D , 

^ -^C2„(r)f^) ^^n{T,^L). (4.22) 



n=0 



Note that there are only even powers of fji/T because of the reflection symmetry, 
Eq. (|4.3p . The leading Taylor coefficient is just the pressure at zero density, and all 
higher coefficients are derivatives evaluated at // = 0, 



Co(T) = ;^(r,M-0), C2n{T) 



" ' {2n)\ 9(f)2" 



(4.23) 

M=0 



These can of course be calculated by standard Monte Carlo techniques. Once the 
coefficients are available, the series for other thermodynamic quantities follow, like the 
quark number density or the quark number susceptibility, respectively, 

n dfl „ n , / u \ 3 



T a(f) T vr. 

H - Hi -2'=" + 12<=.(f)' + 300,(1)' + ^^^ (4,24) 

Clearly, this strategy can be applied to any observable of interest and is well-defined 
computationally. As one needs to calculate coefficient by coefficient, the full functions 
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are only well approximated as long as the series converges sufficiently rapidly, i.e. for 
sufficiently small fi/T. 

In practice, what needs to be evaluated during a simulation are the derivatives of 
a particular observable with respect to chemical potential, 

Using dot Af = cxpTrlnA/, this is converted into expressions like 

aindctA/ / ,(9M 

= Tr Ar^ 



a^lndetA/ _ ^ (^A/^i^^*^^ T f M'^^^Wr^^^^^ 
\ d^i"^ J \ dfj. J 

etc. (4.26) 

The advantage is that the derivatives of the fermion determinant are represented by 
local operators, which can be evaluated using random noise vectors. But it is clear that 
higher order derivatives quickly turn into very complex expressions involving many 
cancellations, and hence are very cumbersome to evaluate with controlled accuracy. 
For a discussion of this point, see ( |Gavai and Gupta, 2005[ ). 

4.5 QCD at imaginary chemical potential 

The hermiticity relation Eq. (|4.12p tells us that the QCD fermion determinant is real 
for imaginary chemical potential ^ = i^i. For such a parameter choice we can simulate 
without sign problem with no more complications than at /x = 0, and moreover we 
can do so with an ensemble that is actually sensitive to chemical potential. In order 
to get back to real chemical potential, we can use once more the Taylor expansion, 

{0){^,)=Y^cui^) , (4.27) 



k=l 



with only even powers for observables without explicit ^-dependence. Note that in 
this case the Monte Carlo results represent the left side of the equation and contain 
no approximation beyond the usual finite size and cut-off effects. Thus, if there are 
sufficiently many and accurate data points, one can test whether the observable is 
well-represented by a truncated Taylor expansion, and if this is the case one may 
analytically continue, fii — > —ifJ-i- 

In order to apply this technique, it is necessary to discuss a few general properties 
of the partition function. From the form of the grand canonical density operator, 
Z = Tr exp — (iJ — ^Q), it is clear that the partition function is going to be periodic 
for imaginary fi. Is this a sensible theory? Modifying the discussion around Eq. (|4.6p . 
imaginary chemical potential is formally equivalent to a real external U{1) field, so 
this parameter choice is well-defined. Next, let us examine the periodicity ( Roberge 
and Weiss, 1986). Since we always have a compactified temporal lattice direction, we 
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can eliminate the quark number term from the aetion and absorb it by the modified 
boundary condition, 

Z^^\ifi^) ^ J DUdctM{0)e-^^, b.c: V(t + iVr, x) = -e'T-?/'(T, x) . (4.28) 

Now let us consider gauge transformed fields 4'^ {x) , (x) , using the large gauge 
transformations discussed in Section [2.5[ i.e. 

5'(r + 7V„a;) =e-'^g'(r,x) . (4.29) 

The measure and the QCD action for /i = are invariant under such a transformation, 
so the new form of the partition function is 

Z^^'>{ifii) = J DU dot M{0)e-'^\ b.c: iA(r + iV^, x) = -e^'^e'^- V(-r, x). (4.30) 

We now observe that 

Since we have obtained one partition function from the other by a gauge transforma- 
tion, the two give completely equivalent descriptions of physics, so we can conclude 
for the QCD partition function 

<'f + '^)^-('f)^ a32, 

Comparing Eqs. (|4.28l4.29p . we see that for discrete values Hi/T — 2TTn/N an imagi- 
nary chemical potential is equivalent to a centre transformation, and Eq. (|4.32p tells 
us that the partition function is symmetric under such transformations, and therefore 
periodic. We have thus established two remarkable properties: the period of the QCD 
partition function in the presence of an imaginary chemcial potential is 2tt/N, and the 
Z{N) centre symmetry is a good symmetry even in the presence of finite mass quarks! 

Let us now specialise to the physical case, 3. We recall from Sec. 12.51 that 

the Polyakov loop closes through the temporal boundary and hence picks up phases 
when changing the centre sector, Eq. (|2.36p . Therefore the phase of the Polyakov 
loop is an observable to identify the different Z'(3)-sectors as the imaginary chemical 
potential is increased. There are Zii) transitions between neighbouring centre sectors 
for all {fj,i/T)c = ^ (n -|- i) ,n = 0, ±1, ±2, .... It has been numerically verified that 
these transitions are first order for high temperatures and a smooth crossover for low 
temperatures ( |de Forcrand and Philipsen, 2002[|D'Elia and Lombardo, 2003[ ). 

Fig. 14.31 (left) shows a schematic phase diagram. The vertical lines represent the 
high temperature, first order Z(3)-transitions. On the other hand, at ^ = we also have 
the chiral/deconfinement transition, which analytically continues to imaginary chem- 
ical potential, as indicated by the dotted lines. As discussed before, its order depends 
on the quark mass combinations. It joins the Z{5) transitions in their endpoint, whose 
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Fig. 4.3 Left: Periodic phase diagram for imaginary chemical potential. Vertical lines cor- 
respond to first order transitions between the different Z(3)-sectors, the arrows indicate the 
corresponding phase of the Polyakov loop. The /i = chiral/deconfinement transition contin- 
ues to imaginary chemical potential, its order depends on Nf and the quark masses. Right: 
Nature of the junction at fixed fii/T — n/3 mod. 2n/3. Solid lines are lines of triple points 
ending in tri-critical points, connected by a Z{2)-lme. 

nature has recently been investigated by explicit simulations for Nf = 2,3 ( D'Elia 
and Sanfilippo, 2009; de Forcrand and Philipsen, 2010[ ). For small and large quark 
masses, where there is a first order chiral and dcconfinement transition, respectively, 
three first order lines join up at this point, rendering it a triple point. For intermediate 
quark masses the quark/hadron transition is only a crossover, and the Z{3) transition 
features an endpoint with 3d Ising universality. Tuning the quark masses interpo- 
lates between these situations, with tri-critical points marking the changes, Fig. 14.31 
(right). If we want to learn about certain observables for real ^ by means of analytic 
continuation, the radius of convergence is given by the closest non-analyticity due to 
phase transitions. At high temperatures and for the chiral/deconfinement transition 
line itself, we are thus limited to the first sector, or /i/T < tt/3. 



4.6 The canonical partition function 

The grand canonical partition function for imaginary chemical potential appears again 
in the formulation of finite density QCD by means of the canonical partition function 
( [Roberge and Weiss, 1986] ) . Consider the Fourier transform of Z(i/ii/r), 




(4.33) 
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In the second line we have made use of the periodicity of the grand canonical partition 
function, Eq. and in the third we have shifted m/T in/T + 2tt/N. The first 

and the third line can only be equal with Z ^ ii 

^=n, n = 0,l,... (4.34) 

Inserting the grand canonical partition function with the quark number operator Q, 
we can write 

= Tr(c-T'5{Q-Q)) , (4.35) 



where we have used the integral representation of the delta-function, 

^(Q-g) = ^/<f)exp(.^^ifc^) . (4.36) 

We thus arrive at the conclusion that Z{Q) corresponds to the QCD partition func- 
tion evaluated with a delta-constraint fixing quark number to be Q. This is just the 
canonical partition function for QCD, 

ZiQ) = ZciQ) . (4.37) 

Since quark number is constrained to be Q = nN, baryon number B = NQ must be 
an integer. 

The descriptions in terms of the canonical or the grand canonical partition func- 
tion are identical in the thermodynamic limit. The grand canonical ensemble can be 
obtained from the canonical one by the fugacity expansion. 



Z(y,T,M) = 5]e^Zc(F,r,B) 

B=l 

/oo 
dpc^^''*Zc(r,p) 
oo 



= lim / dpe-^'^^'^P^-f'P^ , (4.38) 

where B = pV and / has now been defined by the canonical ensemble. It is then 
possible to convert back from baryon number to chemical potential by 

MP) im 

Since the grand canonical partition function at imaginary chemical potential can be 
evaluated by Monte Carlo methods, this offers another approach to deal with finite den- 
sity QCD OAlford, Kapustin and Wilczek, 1999[ ). Numerical data for the grand canon- 
ical partition function can be Fourier transformed numerically to give the canonical 
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partition function at fixed baryon number. However, in this approach the sign problem 
re-enters at the stage of the Fourier transform. The latter clearly has an oscillatory 
integrand and thus the same problem we encountered before. The integration can only 
be done for sufficiently small quark numbers Q, whereas the thermodynamic limit re- 
quires Q ^ oo with a fixed quark number density Q/V = const. Nevertheless, this is 
an interesting alternative approach which does not require reweighting or truncation 
of Taylor series, and thus is explored numerically ( [Kratochvila and de Forcrand, 2006{ 
Alexandru, Faber, Horvath and Liu, 2005| ). 

4.7 Plasma properties at finite density 

Having developed computational tools for finite density, one can apply them to the 
studies discussed in the previous sections and see how finite baryon densities affect 
the screening masses, the equation of state or the static potential. In all those cases 
the influence of the chemical potential is found to be rather weak. These calculations 
appear to be well under control, and we will not further discuss them here. Instead we 
outline recent attempts to determine the QCD phase diagram at finite density, where 
the order of the phase transition is expected to change as fi is increased. Fig. 13.41 
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Towards the QCD phase diagram 



5.1 The critical temperature at finite density 

As in the case of zero density, let us first discuss the phase boundary, Tc{^), before 
deahng with the order of the phase transition. The (pseudo-)critical hne has been 
calculated for a variety of flavours and quark masses using different methods. Its 
computation is most straightforward by reweighting, where one directly evaluates sus- 
ceptibilities of observables, Eq. (|3.33|) . and locates their maxima to identify the critical 
couplings. With the Taylor expansion method, one computes instead the coefficients 
of such susceptibilites, cf. Eq. (|4.24|) . which then is known to some order in ^/T. 
Similarly one can evaluate susceptibilities at imaginary chemical potentials and locate 
their maxima as a function of /i^. 

On the other hand, the (pseudo-) critical gauge coupling can itself be expressed as 
a Taylor series. It was defined as an implicit function to be the coupling for which a 
generalised susceptibility peaks. The implicit function theorem then guarantees that, 
if xiPi m) is an analytic function, so is I3c{li), and hence 



Using some form of the rcnormalisation group beta-function, this can be converted 
into an expansion for the (pseudo-)critical temperature. 



Thus, the Taylor expansion method as well as simulations at imaginary chemical po- 
tential followed by fits to polynomials allow for independent determinations of the 
coefficients, providing valuable cross checks to control the systematics. For a quanti- 
tative comparison one needs data at one fixed parameter set and also eliminate the 
uncertainties of setting the scale. Such a comparison is shown for the critical cou- 
pling in Fig. 15.11 for A^^ = 4 staggered quarks with the same action and quark mass 
m/T K, 0.2. (For that quark mass the transition is first order along the entire curve). 
One observes quantitative agreement up to /x/T w 1.3, after which the different re- 
sults start to scatter. Thus we conclude that all methods discussed here appear to be 
reliable for i^/T < 1. 

For quark mass values close to the physical ones, one observes that the critical 
temperature is decreasing only very slowly with ji on coarse lattices. These calculations 
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Fig. 5.1 Comparison of different methods to compute the critical couplings (Kratochvila 
and de Forcrand, 2006). 

can be repeated on finer lattices. While the cost for reweighting becomes formidable, 
the other methods should allow for a reliable determination of the leading coefficients 
in the continuum limit, and thus the physical phase boundary between the hadron and 
quark gluon plasma phase for /i < T, in the near future. 



5.2 The QCD phase diagram for /i 7^ and the critical point 

As in the case of /i = 0, a determination of the order of the transition, and hence 
the search for the critical endpoint, is much harder, and we begin by discussing the 
qualitative picture. If a chemical potential is switched on for the light quarks, there 
is an additional parameter requiring an additional axis for our phase diagram char- 
acterising the order of the transition. Fig. 13.81 This is shown in Fig. 15.21 where the 
horizontal plane is spanned by the n = phase diagram in m^, m„_d and the vertical 
axis represents fi. The critical line separating the first order section from the crossover 
will now extend to finite /x and span a surface. A priori the shape of this surface is not 
known. However, the expected QCD phase diagram corresponds to the left scenario 
of Fig. 15.21 The first order region expands as /i is turned on, so that the physical 
point, initially in the crossover region, eventually belongs to the chiral critical surface. 
At that chemical potential fiE, the transition is second order: that is the QCD chiral 
critical point. Increasing fi further makes the transition first order. In this scenario the 
transition is generally strenghtened with real chemical potential. A completely differ- 
ent scenario arises if instead the transition weakens and the first order region shrinks 
as /X is turned on. In that case the physical point remains in the crossover region and 
there is no chiral critical point for moderate /x. Fig. 15.21 (right). 

There are then two different strategies to learn about the QCD phase diagram. One 
can fix a particular set of quark masses and for that theory switch on and increase 
the chemical potential to see whether a critical surface is crossed or not. Sec. 15.31 
Alternatively, Sec. 15. 41 discusses how to start from the known critical line at /x = and 
study its evolution with a finite /Lt. 
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Real world Real world 

Heavy quarks n Heavy quarks 




Fig. 5.2 Upper panel: The chiral critical surface in the case of positive (left) and negative 
(right) curvature. If the physical point is in the crossover region for = 0, a finite n chiral 
critical point will only arise in the scenario (left) with positive curvature, where the first-order 
region expands with /i. Note that for heavy quarks, the first-order region shrinks with /i, 
cf. Sec. 15.51 Lower panel: phase diagrams for fixed quark mass (here Nf = 3) corresponding 
to the two scenarios depicted above. 

5.3 Critical point for fixed masses: reweighting and Taylor 
expansion 

Reweighting methods at physical quark masses get a signal for a critical point at 
fi^ ~ 360 MeV dFodor and Katz, 2004| . In this work x A lattices with L = 6 - 12 
were used, working with the standard staggered fermion action. Quark masses were 
tuned to rriu^d/Tc ~ 0.037, TOj/Tc ss 1, corresponding to the mass ratios mT^/mp « 
0.19, mTr/rriK ~ 0.27, which are close to their physical values. A Lee- Yang zero analysis 
was employed in order to find the change from crossover behaviour at /x = to a first 
order transition for /i > /i^;. This is shown in Fig. 15.31 For a crossover the partition 
function has zeroes only off the real axis, whereas for a phase transition the zero moves 
to the real axis when extrapolated to infinite volume. For a critical discussion of the 
use of Lee- Yang zeros in combination with reweighting, see ( |Ejiri, 2006[ ). A caveat of 
this calculation is the observation that the critical point is found in the immediate 
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Fig. 5.3 Imaginary part of tlie Lee- Yang zero closest to tiie real axis as a function of 
cliemical potential ( |Fodor and Katz, 2004p . 

neighbourhood of the onset of pion condensation in the phase quenched theory, which 
is w here sign problem becomes maximaUy severe (fSphttorff, 2005[ Han and Stephanov, 
2008). cf. the discussion in Sec. 14.31 Therefore, one would hke to confirm this result 
with independent methods. 

In principle the determination of a critical point is also possible via the Taylor 
expansion. In this case true phase transitions will be signalled by an emerging non- 
analyticity, or a finite radius of convergence for the pressure series about fi = 0, 
Eq. (j4.22p . as the volume is increased. The radius of convergence of a power series 
gives the distance between the expansion point and the nearest singularity, and may 
be extracted from the high order behaviour of the scries. Possible definitions are 



p, r = lim pn,r„ 



with 



Co 

C2n 



l/2n 



C2n 



C2n+2 



1/2 



(5.3) 



General theorems ensure that if the limit exists and asymptotically all coefficients 
of the series are positive, then there is a singularity on the real axis. More details as 
well as previous applications to strong coupling expansions in various spin models can 
be found in ( [Guttman, 1989[ ). In the series for the pressure such a singularity would 
correspond to the critical point in the (/i, r)-plane. The current best attempt is based 
on four consecutive coefficients, i.e. knowledge of the pressure to eighth order, and a 
critical endpoint for the Nf = 2 theory was reported in (Gavai and Gupta, 2005). 
There are also difficulties in this approach. Firstly, there are different definitions for 
the radius of convergence, which are only unique in the asymptotic limit, but differ 
by numerical factors at finite n. Furthermore, the estimated Pn,rn at a given order 
is neither an upper nor a lower bound on an actual radius of convergence. Finally, at 
finite orders the existence of a finite radius of convergence is a necessary, but not a 
sufficient condition for the existence of a critical point. For example, one also obtains 
finite estimates for a radius of convergence from the Taylor coefficients of the hadron 
resonance gas model, even though that model does not feature a non-analytic phase 
transition. 
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5.4 The change of the chiral critical line with ^ 

Rather than fixing one set of masses and considering the effects of /x, one may map 
out the critical surface in Fig. I5.2l bv measuring how the ^ = critical boundary line 
changes under the influence of ^. For example, for a given value of the strange quark 
mass the corresponding critical m^^d has again a Taylor expansion, 

+ , (5.4) 



The curvature of the critical surface in lattice units is directly related to the behaviour 
of the Binder cumulant via the chain rule. 



darric dB^ ( dB^ 



d{a^Y d{a^Y \dam 



(5.5) 



and similar expressions for the higher derivatives. While the second factor is sizeable 
and easy to evaluate, the /.t-dependence of the cumulant is excessively weak and re- 
quires enormous statistics to extract. In order to guard against systematic errors, this 
derivative can be evaluated in two independent ways. One is to fit the corresponding 
Taylor series of ^4 in powers of ji/T to data generated at imaginary chemical po- 
tential, the other to compute the derivative directly and without fitting via the finite 
difference quotient, 

dB, _ ^.^ B,ia,) Bm_ (5.6) 



d{a^Y (a^i)2-i.o (a/i)2 

Because the required shift in the couplings is very small, it is adequate and safe to use 
the original Monte Carlo ensemble for am'^(O), /i = and reweight the results by the 
standard Ferrenberg-Swendsen method. Moreover, by reweighting to imaginary /i the 
reweighting factors remain real positive and close to 1. 

On coarse Nr ~ 4 lattices, the first two coefficients in Eq. (|5.4[) are found to be 



negative, ( de Forcrand and Philipsen, 2008 1 , hence the region featuring a first order 



chiral transition is shrinking when a real chemical potential is turned on. This implies 
that, for moderate /i<T, there is no critical point belonging to the chiral critical 
surface. An open questions here is what the higher order corrections are. For example, 
a sudden change of behaviour as in Fig. 15.31 would be difficult to capture with only 
few Taylor coefficients. Also, there might be a chiral critical point at larger chemical 
potentials and finally we do not yet know yet how the curvature of the critical surface 
behaves in the continuum limit. In addition, these findings do not exclude a critical 
point that is not associated with the chiral critical surface. 

5.5 The change of the deconfinement critical line with /i 

In light of these results, it is interesting to compare with the situation in the heavy 
quark corner of the schematic phase diagram. Fig. 15.21 As the quark mass goes to 
infinity, quarks can be integrated out and QCD reduces to a gauge theory of Polyakov 
lines. At a second order phase transition, universality allows us to neglect the details 
of gauge degrees of freedom, so the theory reduces to the 3d three-state Potts model. 
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Fig. 5.4 The critical heavy quark mass separating first order from crossover as a function 
of fi^ from the Potts model ( |Kim, de Forcrand, Kratochvila and Takaishi, 2006[ |. 

which is in the appropriate 3d Ising universahty class. Hence, studying the three-state 
Potts model should teach us about the behaviour of QCD in the neighbourhood of 
the deconfinement critical line separating the quenched first order region from the 
crossover region. 

Here we are interested in simulations at small chemical potential. In this case, the 
sign problem is mild enough for brute force simulations at real /i to be feasible. In 



(Kim, de Forcrand, Kratochvila and Takaishi, 2006), the change of the critical heavy 
quark mass is determined as a function of real as well as imaginary fi, as shown in 
Fig. 15.41 Note that M'^{^) rises with real chemical potential, such that the first order 
region in Fig. 13.81 shrinks as finite baryon density is switched on. The calculation also 
gives some insight in the problem of analytic continuation: Fig. 15.41 clearly endorses 
the approach with Mdp-) passing analytically through = 0. However, high order fits 
might be required in practice in order to reproduce the data on both sides of = 0. 

The same qualitative behaviour is observed in a combined strong coupling and 
hopping parameter expansion ( [Langelage and Philipsen, 2010a[ ). While such a calcu- 
lation is only accurate for coarse lattices, the universal features of the critical surface 
should be cut-off independent. Thus, QCD with heavy quarks is an example of the 
non-standard scenario discussed in the previous section. 



5.6 The critical surfaces at imaginary fi 

The chiral and deconfinement critical surfaces, upper Fig. 15.21 (right), continue to 
imaginary fx, which has been used to determine their curvature. Rather than studying 
the neighbourhood of = with an aim to analytic continuation, it is also illuminating 
to follow those surfaces into the imaginary regime, until the phase boundary to the 
next Z(3)-sector, Fig. 14.31 (left'), is hit. This has already been done in the case of the 
deconfinement transition within the Potts model, cf. the curve at negative Fig. 15.41 
The meeting point of the deconfinement critical line with the critical endpoint of the 
Z(3) transition is then tri-critical, and corresponds precisely to the tri-critical point 
in the large mass region of Fig. 13.81 (right). 
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Fig. 5.5 Left: Nature of the Z(3)-transition endpoint at ^i/T — ivr/S. (de Forcrand and 
Philipsen, 2010). Right: Deconfinement critical surface for real and imaginary fi. 

For three flavours of quarks and varying quark masses, one can now draw a 
schematic diagram analogous to Fig. 13. 4[ showing the nature of the junction at /i = 
i7rT/3. At the time of writing, only the two tri-critical points on the Nf = 3 diagonal 
and the light mass one for Nf = 2 have been determined on coarse lattices. A natural 
extension would then be that there are two tri-critical lines which delimit the quark 
mass sections for which the junction corresponds to triple points. The area in between 
corresponds to a second order endpoints of the Z(3)-transition, with no connection to 
the chiral/deconfinement transition which is merely a crossover there. The tri-critical 
lines represent the boundaries in which the chiral and deconfinement critical surfaces 
end at imaginary fi ~ i-kT /i, before they periodically repeat. On coarse lattices they 
are found "outside" the corresponding /i = critical lines, which demonstrates an 
increase of the first order area in the imaginary direction. In fact, for heavy quarks 
this change is monotonic. The functional form of the curve in Fig. 15. 4[ and hence the 
surface in Fig. 15.51 (right), was found to be determined by tri-critical scaling for the 
whole range up to real fi ^ T. This would explain the curvature of the critical sur- 
faces, i.e. the weakening of the chiral and deconfinement transitions with real chemical 
potential. 



5.7 Discussion 

It thus appears that first order transitions are weakened when a chemical potential 
for baryon number is switched on. The same observation is made for finite isospin 



chemical potential (Kogut and Sinclair, 2008). Note however, that most of these com 



putations are done on Nj- = 4, 6 lattices, where cut-off effects appear to be larger than 
finite density effects. Hence, definite conclusions for continuum physics cannot yet be 
drawn. A general finding is the steepness of the critical surface, making the location 
of a possible critical endpoint extremely quark mass sensitive, and hence difficult to 
determine accurately. Furthermore, the chiral first order region is also observed to 
shrink with decreasing lattice spacing, such that in Fig. 13.81 mj"'' < mf''''*. In this 
case the chiral critical surface we have been discussing here might not be responsible 
for a possible critical point, regardless of its curvature, but another surface emanating 
from the putative 0(4)-chiral limit, or one connected to finite density physics, or 
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Possibilities and uncertainties abound! To weed them out we need to conclusively un- 
derstand the situation in the Nf ~ 2 chiral limit. Finally, even if those issues get 
settled, all the methods discussed here are only valid for < T. Thus, learning about 
high density QCD, relevant for compact stars, requires entirely new theoretical meth- 
ods. In summary, it remains a challenging but also most exciting task to settle even 
the qualitative features of the QCD phase diagram. 
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